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Iterative  solution  of  multiple  right  hand  side  matrix  equations  for  frequency  domain 

monostatic  radar  cross  section  calculations 

M  D  Pocock,  S  P  Walker 

Computational  Mechanics,  Imperial  College  of  Science  Technology  and  Medicine,  London  SW7  2BX 

s.p.walker@ic.ac.uk 


Abstract 

Monostatic  res  characterisation  using  integral 
equation  methods  in  the  frequency  domain  requires 
the  solution  of  very  large  matrix  equations  with 
multiple  right  hand  sides.  Although  costly  for  a 
single  right  hand  side,  direct  methods  are 
attractive  in  that  subsequent  right  hand  sides  are 
very  cheap.  Iterative  methods  are  much  cheaper 
for  a  single  right  hand  side,  but  if  the  whole 
solution  must  be  repeated  for  each,  they  become 
much  more  expensive.  We  investigate  here  the 
performance  of  a  simple  modification  to  the  GCR 
algorithm,  which  allows  solutions  for  an 
essentially  unlimited  number  of  right  hand  sides  to 
be  obtained  for  a  modest  multiple  of  the  cost  of  the 
first.  For  the  cases  investigated,  with  up  to  360 
right  hand  sides  on  bodies  up  to  15  wavelengths 
long,  with  matrices  up  to  20,440  by  20,440  in  size, 
this  multiple  was  below  -10.  Costs  seem  to  rise 
with  the  number  of  right  hand  sides  till  the  surface 
field  is  in  some  sense  characterised,  and  thereafter 
subsequent  illumination  angles  are  essentially  free. 
An  investigation  of  cost  scaling  on  a  set  of  spheres, 
ranging  from  -1  to  7  wavelengths  in  diameter,  seems 
to  indicate  the  cost  of  full  monostatic 
characterisation  to  scale  with  about  the  fifth 
power  of  frequency. 

1.  Introduction 

There  are  two  main  cost  components  in  the  solution 
of  large  scattering  problems,  such  as  res  evaluation, 
via  frequency  domain  integral  equation  methods. 
The  first1,  scaling  with  frequency/  to  the  fourth 
power,  is  the  cost  of  formation  of  the  dense  matrix 
describing  interactions  between  parts  of  the 
scatterer.  The  second  is  the  cost  of  solving  the 
resulting  matrix  equation. 

For  the  direct  solution  schemes  most  usually 
employed,  this  solution  cost  scales  with  N3,  where 
N  is  the  order  of  the  matrix  equation,  or 
equivalently  with/6.  For  all  but  small  problems 
matrix  solution  is  thus  the  dominant  cost.  Because  of 
this  large  cost,  there  has  been  increasing  attention 
paid  of  late  to  the  use  of  iterative  methods  for 


matrix  solution2'8.  It  is  generally  concluded  from 
such  studies  that  iterative  schemes  require  much 
less  time  to  compute,  often  by  an  order  of  magnitude 
or  more,  than  would  be  taken  to  solve  the  same 
system  by  direct  methods.  Cost  reduction  by  a  large 
multiple  is  obviously  welcome;  even  more  so  would 
be  if  that  multiple  itself  increased  with  problem 
size,  corresponding  to  a  reduction  in  cost  scaling 
below  the  /6  of  the  direct  approach.  This  has 
indeed  been  suggested7,  but  the  evidence  is  as  yet 
inconclusive8. 

For  practical  radar  cross  section  analyses,  it  is 
necessary  to  determine  the  scattered  field  caused  by 
fields  incident  from  a  (large)  number  of  directions. 
Direct  solution  methods  are  then  attractive.  Once 
the/6  cost  of  matrix  factorisation  has  been  paid, 
solutions  for  large  numbers  of  incident  waves  (right 
hand  side  vectors  in  matrix  equation  terms)  can 
each  be  found  at  a  cost  scaling  with/4.  Even  if  the 
required  number  of  right  hand  sides  scales  with  the 
second  power  of  frequency,  the  overall  cost  of  a  full 
monostatic  evaluation  still  only(!)  scales  with 
frequency  to  the  sixth  power. 

This  interest  in  multiple  illumination  angles  is  a 
discouragement  to  the  use  of  iterative  methods.  One 
approach  using  them  is  simply  to  start  the  solution 
afresh  for  each  right  hand  side.  Naturally,  if  more 
than  a  few  right  hand  sides  are  required,  any  cost 
saving  over  the  direct  approach  is  lost.  If  the 
required  number  of  right  hand  sides  scales  with  the 
second  power  of  frequency,  the  overall  cost  then 
scales  with  something  approaching/8. 

We  present  here  methods  which  seem  able  to 
provide  iterative  solutions  for  additional 
monostatic  illumination  angles  (right  hand  sides) 
at  essentially  no  cost,  once  sufficient  right  hand 
sides  have  been  analysed  to  characterise,  in  some 
sense,  the  surface  field. 

There  have  been  a  few  attempts  to  address  this 
multiple  right  hand  side  issue  for  iterative 
solutions,  which  we  will  discuss  prior  to  describing 
the  present  approach. 
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In  general  the  quality  of  the  initial  guess  has  only  a 
modest  affect  on  the  number  of  iterations  required 
for  convergence  of  an  iterative  method2-4,7. 
Consequently,  the  approach  of  using  the  solution 
from  one  incident  wave  as  the  starting  point  for  the 
next  offers  little  benefit. 

For  the  sparse  matrices  resulting  from  a  finite 
element  discretisation.  Smith  and  colleagues6 
effectively  expanded  several  solutions  in  terms  of  a 
single  set  of  search  vectors,  with  some  success. 

In  a  very  recent  paper9  Boyse  and  Seidl  use  the 
GMRES  algorithm  (and  see  references  cited  in  Boyse 
for  more  details  of  the  block  GMRES  and  the 
multiple  right  hand  side  variant  thereof).  In 
essence  the  problem  is  first  solved  simultaneously 
for  a  modest  number  of  right  hand  sides,  distributed 
uniformly  over  the  span  of  right  hand  sides  of 
eventual  interest,  using  a  block  GMRES  approach. 
Whilst  this  costs  more  than  for  a  single  right  hand 
side,  it  is  not  a  large  multiple  of  the  cost. 
Intermediate  values  are  then  found,  using  the 
orthonormal  basis  for  the  Krylov  sub-space  which 
was  found  during  the  solution  for  the  initial  right 
hand  sides.  It  was  found  necessary  to  ensure  that 
the  number  of  right  hand  sides  solved  for  initially 
was  carefully  chosen.  This  number  of  right  hand 
sides  was  noted  as  needing  to  be  sufficient  to 
represent  in  some  sense  the  RCS  distribution  sought. 
An  angular  separation  such  as  to  allow  rather  more 
than  two  samples  per  wavelength  was  suggested.  If 
this  is  not  observed,  intermediate  solutions  found 
subsequently  tended  to  have  significant  errors.  Too 
many  initial  right  hand  sides,  however,  caused 
slow  convergence  of  the  block  GMRES  solution. 
Overall,  significant  reduction  in  time  relative  to  a 
direct  solution  was  observed. 

The  GCR  (generalised  conjugate  residual) 
algorithm10  was  used  by  Soudais5  for  solving  the 
matrix  equations  resulting  from  analysis  of 
scattering  from  a  -1  wavelength  mixed  dielectric- 
perfect  conductor  target.  A  finite  element  and 
symmetric  Stratton-Chu  treatment  were  combined, 
giving  a  matrix  system  which  was  symmetric,  and 
sparse  in  the  finite  element  regions.  Multiple  right 
hand  sides  were  tackled  by  what  seemed  to  be  a 
very  effective,  and  attractively  simple, 
modification  to  the  normal  GCR  algorithm.  As 
with  most  iterative  methods,  the  solution  is 
changed  at  each  iteration  by  moving  a  particular 
distance  (the  step  length)  in  some  direction  (the 
search  direction).  The  computationally  expensive 
part  of  the  algorithm  is  the  finding  of  this  search 


direction.  In  the  modified  algorithm,  the  search 
direction  is  found  for  only  that  right  hand  side 
currently  exhibiting  the  largest  residual.  The 
solutions  to  all  right  hand  sides  are  advanced  in 
this  direction,  but  for  a  different  (and  optimal)  step 
length  in  that  direction  for  each. 

In  the  sections  which  follow  we  will  extend  the 
application  of  this  approach,  to  the  dense  and 
unsymmetric  matrices  which  result  from  a  normal 
integral  equation  discretisation,  and  to  analysing 
scattering  from  multi-wavelength  'stealthy' 
targets.  A  major  aspect  of  interest  will  be  the 
interaction  of  the  number  of  iterations  required  (a 
measure  of  the  computational  work),  with  the 
number  of  right  hand  sides  analysed,  and  the  body 
size  in  wavelengths.  It  is  this  latter  which 
determines  the  'jaggedness*  with  angle  of  the 
monostatic  radar  cross  section,  and  consequently  the 
number  of  illumination  angles  required  to 
characterise  it.  If  the  number  of  iterations  continues 
to  rise  in  proportion  to  the  number  of  illumination 
angles,  till  the  point  that  the  response  is  fully 
characterised,  there  may  be  no  reduction  in  cost 
scaling;  if  otherwise,  there  could  be. 

In  the  next  section  we  will  summarise  the  integral 
equation  formulation,  and  outline  the  normal  GCR 
algorithm,  and  the  modification  to  handle 
multiple  right  hand  sides.  Section  three  will 
present  results  from  its  application  to  a  number  of 
res  problems. 

2.  Formulation 

The  formulation  is  well  known,  and  only  a  brief 
description  will  be  given  here.  The  MFIE  for 
scattering  from  a  smooth,  closed  perfect  conductor  is 

}H(r)  =  ff"c(r)  + 

1  „  f  1  (!) 

oQ. 

where  the  integrations  to  obtain  the  field  at  surface 
location  r  are  over  the  rest  of  the  surface  of  the 
scatterers'  andr,  withn  the  unit  normal,  and  R  the 
vector r-r1 .  The  incident  wave  is  H'^.  We  employ  a 
curvilinear,  isoparametric  discretisation,  with 
Gaussian  quadrature.  Fuller  details  are  given 
elsewhere11;  the  result  is  a  (complex)  matrix 
equation  of  the  form: 

[A]H  =  H  "c  (2) 
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This  is  to  be  solved  with  tens  or  possibly  hundreds 
of  incident  wave  vectors  H^. 

We  use  the  unpreconditioned  GCR  algorithm, 
which  we  first  outline  in  its  normal  ‘single  right 
hand  side'  form. 


To  solve 
Ax  =  b 
Initialise: 

x(0)=0 

r<°>=-b 

p(0)=_r(°) 

For  each  iteration  k,  fork  =  0, 
Calculate  the  step  length  a: 


a 


(k)  _  . 


,Ap«) 


Ap<‘> 

Increment  the  solution  vector 
x(t+1)=xw+a«pW 
and  residual  vector: 
r<‘«)=r(*>+aWApM 

Check  if  converged  sufficiently:  stop  if 

,(*+ 1)|| 

<10~m 


M 


(3) 


(4) 


(5) 


(6) 


(7) 


(8) 


withm  typically  in  the  range  2  to  8. 

I  f  not,  calculate  search  direction  vector  and  matrix 
vector  product  for  the  next  iteration: 


<7 


(*+i) 


(Ar(*+1),Ap(/)) 


Ap 


.(0 


fori  = 


(9) 


p(*+i)  =  _r(*+i) + Yof+1)p(,) 


!=1 


(10) 


Ap(fc+1)  =  -Ar(fc+1>  +Jja(ik+1)Ap^ 


i=l 


(11) 


As  well  as  the  main  system  matrix,  we  see  that  the 
sets  of  vectors  p(k)  and  Ap(k)  must  be  stored. 

The  modifications  required  to  solve  for  multiple 
right  hand  sides  are  very  simple.  In  the 
initialisation  stage  a  residual  and  initial  search 
direction  is  obtained  for  each  right  hand  side. 

At  equation  (5)  a  step  length  is  calculated  for  each 
right  hand  side,  and  the  solution  and  residual 
vectors  incremented.  Between  (7)  and  (8)  the  largest 
(normalised)  residual  is  found,  by  searching 
through  the  residuals  associated  with  all  right 
hand  sides. 

Assuming  it  fails  the  test  (8),  a  single  search 
direction  is  calculated  in  (9)  and  (10),  to  suit  that 
right  hand  side  found  to  have  the  largest  residual. 
This  step,  involving  matrix  vector  multiplication, 
is  the  one  where  the  principal  computation  lies.  A 
different  step  length  for  each  right  hand  side  is 
then  calculated  in  (5),  which  is  used  with  the 
common  search  direction  to  increment  each  solution 
in  (6).  We  see  that  at  each  iteration  the  parameters 
selected  most  suit  that  right  hand  side  currently 
furthest  from  solution.  As  a  result  that  right  hand 
side  improves  rapidly,  with  some  other  right  hand 
then  side  taking  its  place  as  the  current  worst.  Thus 
are  all  right  hand  sides  shepherded  to  solution 
more  or  less  in  step. 

Each  right  hand  side  is  seen  to  increase  by  one  the 
number  of  vector-vector  manipulations  required  a  t 
each  iteration  (multiplications  in  equation  (5)),  and 
vector-vector  additions  in  equations  (6),  (7)  and 
(8)).  Storage  is  increased  by  the  need  to  store  the 
evolving  solution  and  residual  for  each  right  hand 
side.  For  problems  of  practical  interest,  these 
increases  are  a  modest  fraction  of  the  requirements 
for  a  single  right  hand  side.  As  the  solutions  are 
attained,  it  would  be  possible  to  continue  to  iterate 
only  for  those  whose  residual  is  not  yet  below  the 
specified  tolerance.  Although  this  would  save  some 
cost,  the  saving  would  be  small,  and  the  refinement 
has  not  been  implemented. 

3.  Demonstrations 

In  the  sections  which  follow  we  present  results  of 
the  analysis  of  targets  illuminated  by  many 
different  incident  waves.  Of  particular  interest 
will  be  the  variation  of  the  total  number  of 
iterations  with  the  number  of  illumination  angles, 
and,  for  practical  purposes,  the  variation  of  the  res 
over  that  range  of  illumination  angles.  Does  the 
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number  of  iterations  rise  with  illumination  angles 
as  long  as  'new'  information  is  being  gained? 

3.1  Equatorial  scans 

We  consider  first  a  monostatic  equatorial  scan  of 
the  NASA  almond12,  illuminated  at  7  GHz,  making 
it  ~6  wavelengths  long.  Discretisation  employed  an 
average  nodal  separation  of  about  1/9  of  a 
wavelength,  and  throughout  a  termination  residual 
of  lCT4  was  employed.  This  is  a  value  recently 
suggested13  as  appropriate  for  this  discretisation, 
for  the  curvilinear  isoparametric  modelling 
employed. 

Figure  1  shows  the  variation  in  the  number  of 
iterations  required  with  the  number  of  illumination 
angles.  The  monostatic  scan  was  computed  as  a 
single  run,  with  respectively  1,  5,  10,  ...  up  to  180 
illumination  angles  (right  hand  sides)  in  turn. 
Cases  were  studied  with  these  illumination  angles 
distributed  uniformly  over  180°,  90° ,  45°  and  15°, 
forming  the  lines  shown  on  the  figure. 

We  see  that  in  each  case  the  number  of  iterations 
required  initially  rises  with  the  number  of 
illumination  angles  examined.  For  the  180^  case, 
this  rise  is  essentially  complete  some  time  before 
~30  illumination  angles,  by  when  about  six  times  as 
many  iterations  are  required  as  are  required  for  a 
single  illumination  angle.  By  this  point  we  are 
gaining  ~30  solutions  for  a  multiple  of  ~6  in  the 
work  required.  From  then  on  solutions  for  additional 
angles  are  obtained  without  further  iterations;  the 
maximum  180  considered  are  similarly  gained  with 
this  same  multiple  of  ~6  over  the  cost  of  one.  (As 
discussed  in  section  2  above,  there  is  still  some 
modest  angle-dependent  cost,  in  the  vector 
manipulations  of  (5)  -  (8).  This  is  very  small  in 
practice,  and  will  be  neglected  in  subsequent 
discussion.)  Similar  comments  can  be  made  about 
cases  with  the  illumination  angles  distributed  over 
a  smaller  range.  As  the  extent  is  reduced,  in  each 
case  the  'plateau'  in  iterations  occurs  at  a  lower 
number  of  iterations,  and  after  fewer  illumination 
angles.  For  example,  for  a  15°  degree  range,  the 
plateau  in  number  of  iterations  is  reached  after 
about  10  illumination  angles  compared  to  the  ~30  of 
the  180P  case.  This  is  consistent  with  the  number  of 
iterations  being  a  function  of  the  amount  of 
information  being  sought;  reducing  the  range  reduces 
the  amount. 

There  is  an  interesting  parallel  here  with  the 
observations  of  Boyes  and  Seidl9.  Using  their  very 
different  approach,  they  found  costs  to  rise 


essentially  linearly  with  number  of  illumination 
angles  till  the  response  was  partially 
characterised,  and  that  subsequent  illumination 
angles  were  essentially  free. 

However,  the  response  in  that  reference  was 
considered  primarily  in  terms  of  res,  although  the 
actual  unknown  for  which  we  are  solving  a  set  of 
equations  is  the  surface  field.  The  monostatic  res, 
the  quantity  actually  sought,  and  the  surface  field 
at  any  particular  surface  location,  are  very 
different  functions  of  incident  illumination  angle. 
The  (monostatic)  res  generally  is  a  much  less  smooth 
a  function  of  incident  illumination  angle  than  is  the 
surface  field.  (For  simple  geometries  it  can  of  course 
be  much  more  smooth,  as  the  case  of  the  monostatic 
res  of  a  multi-wavelength  sphere  attests.) 

We  find  some  empirical  evidence  that  the  plateau 
is  reached  when  the  density  of  illumination  angles 
is  such  as  to  represent  reasonably  the  variation  of 
surface  field  with  illumination  angle.  In  figure  2  is 
shown  an  indication  of  the  variation  of  surface 
field  with  illumination  angle  for  the  six 
wavelength  almond,  by  selecting  a  node  (at  -4.1  ,0, 
0)  and  component  (the  real  part  of  the  y  component) 
a  t  random.  The  locations  of  30  uniformly  spaced 
illumination  angles  are  marked  on  this  figure,  and 
it  could  be  argued  that  results  from  only  this  many 
illumination  angles  provide  a  reasonably  good 
representation  of  the  surface  field  variation. 

The  corresponding  res  is  shown  in  figure  3,  with  the 
solid  line  obtained  by  illumination  from  180  angles, 
and  the  30  locations  again  marked.  (This  present 
paper  is  not  concerned  with  characterisation  of  an 
res  code,  but  for  completeness  this  figure  also  shows 
the  measured  res12,  albeit  with  the  measured  values 
obtained  from  an  enlarged  photocopy  of  the 
published  measurements.  As  is  seen,  agreement  of 
the  computed  values  with  this  is  good.)  We  see 
that  much  of  the  detailed  structure  of  the  res  is  not 
revealed  by  using  only  these  30  angles;  lobes  cn  the 
RCS  plot  are  generally  seen  to  be  significantly 
narrower  than  on  that  of  the  surface  field. 
However,  with  the  30  illumination  angles 
corresponding  to  the  start  of  the  plateau,  this 
detail  is  obtained  essentially  free.  (We  have  not 
investigated  the  point,  but  this  might  be  interesting 
to  consider  in  the  context  of  the  'monostatic  - 
bistatic  approximation'14.  Here  relatively  widely 
separated  illumination  angles  can  be  used  to 
provide  a  good  approximation  to  the  finely 
sampled  monostatic  response.  It  may  indeed  be  that 
similar  criteria  apply.) 
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Timings  may  be  of  interest.  On  a  Dec  Alpha  600 
workstation,  the  matrix  took  about  16  minutes  to 
form,  and  25  minutes  to  solve  for  each  illumination, 
giving  a  total  time  of  75  hours  for  a  180  illumination 
angle  characterisation.  Using  the  multiple  right 
hand  side  approach,  these  same  180  solutions  were 
obtained  in  about  3  horns. 

As  a  further  example,  figure  4  (inset)  shows  the 
monostatic  res  of  a  11:1  cylindrical  dipole,  with 
hemispherical  end  caps,  discretised  with  3026 
nodes,  and  illuminated  such  as  to  make  it  -15 
wavelengths  long.  This  figure  was  drawn  from 
results  evaluated  with  360  illumination  angles  0.5° 
apart,  and  shows  approximately  60  distinct  peaks 
over  the  180°.  It  is  not  shown,  but  the  variation  of 
surface  field  with  illumination  angle  is  not 
surprisingly  very  much  smoother  than  the  res  cn 
this  15  wavelength  body;  indeed  it  is  very  smooth 
for  many  surface  locations,  with  -15  peaks  being 
the  maximum  found.  Figure  4  shows  the  variation  of 
the  number  of  iterations  with  the  number  of 
illumination  angles.  Again,  we  see  a  distinct 
plateau  occurring,  here  after  -90  illumination 
angles.  This  is  sufficient  to  represent  reasonably  the 
surface  field  variation,  but  not  to  characterise  the 
equatorial  res  variation.  Here  the  full  equatorial 
monostatic  result,  employing  360  illumination 
angles,  is  obtained  for  -4  times  the  cost  of  a  single 
iterative  solution.  For  this  size  of  matrix,  6052  by 
6052  (complex),  we  generally  find  a  single  iterative 
solution  to  cost  rather  more  than  an  order  of 
magnitude  less  than  a  single  direct  solution. 

3.2  Near  head-on  res 

Whilst  it  is  conventional  to  analyse  equatorial 
sweeps  as  above,  practical  interest  may  probably  be 
concentrated  on  a  relatively  small  solid  angle 
centred  around  head-on,  possibly  biased  towards 
illumination  from  slightly  below.  We  have 
analysed  this  same  6  wavelength  almond,  with 
illumination  in  the  range  0  to  16°  vertically  and  0  to 
16°  horizontally  from  head-on  (where  symmetry 
naturally  makes  only  this  one  quadrant  necessary). 
A  uniform  increment  in  each  angular  coordinate  was 
used,  with  computations  made  with  16  by  16  (256) 
illumination  angles,  8  by  8, 4  by  4,  3  by  3  and  1  (head 
on)  illumination  angle.  The  inset  in  figure  5  shows  a 
representation  of  the  res,  plotted  from  the  256 
illumination  angle  result.  Figure  5  shows  the 
variation  in  the  number  of  iterations  with  the 
number  of  illumination  angles,  again  exhibiting  the 
characteristic  'plateau',  here  at  about  7  times  the 


number  of  iterations  required  for  a  single 
illumination  angle.  ° 

3.3  Cost  scaling 

As  noted  earlier,  there  is  ro  clear  evidence 
regarding  cost  scaling  for  single  illumination  angle 
(single  right  hand  side)  solutions  via  iterative 
methods.  The  position  is  naturally  made  more 
complicated  once  multiple  illumination  angles  are 
included.  We  can  present  here  some  empirical 
evidence,  but  only  tentative  observations  and 
conclusions  can  be  drawn. 

As  discussed  elsewhere8,  fineness  of  discretisation 
(expressed  in  terms  of  degrees  of  freedom  per 
incident  wavelength)  itself  can  influence  the 
number  of  iterations  required.  We  will  here  employ 
spherical  scatterers  of  a  range  of  sizes,  as  uniform 
discretisation  is  difficult  to  ensure  on  (say)  the 
almond. 

Whilst  obviously  it  is  machine  dependant,  some 
actual  times  might  be  helpful.  All  jobs  were  run  on  a 
96Mb  SGI  Indy  R5000.  The  matrix  for  the  biggest 
mesh,  20440  (single  precision  complex)  x  20440  in 
size,  occupied  3.2  Gb.  This  was  formed  once,  and 
read  in  for  each  iteration.  Each  such  iteration  took 
-10  minutes,  of  which  almost  half  was  reading  from 
disk.  At  the  plateau  of  figure  6,  mentioned  below, 
some  680  iterations  were  required,  corresponding  to 
a  time  of  110  horns  for  the  full  characterisation. 

We  plot  in  figure  6  the  number  of  iterations  required 
versus  illumination  angles  for  a  180°  scan  of  a  series 
of  spheres.  These  range  in  size  from  1.6  wavelengths 
in  diameter  (1060  by  1060  matrix)  to  7.12 
wavelengths  (20440  by  20440  matrix).  They 
display  behaviour  qualitatively  identical  to  that 
of  the  almond,  with  a  flat  plateau  in  number  of 
iterations  required  being  reached  after  only  a 
modest  number  of  illumination  angles.  Here, 
though,  the  monostatic  res  is  of  course  characterised 
fully  by  a  single  illumination,  whereas  the 
computed  surface  field  distribution  is  dependant  on 
the  illumination  angle,  lending  further  evidence  to 
the  observation  that  it  is  the  latter  which  is  the 
determinant  of  computational  work. 

In  figure  7  we  plot  the  variation  of  the  number  of 
iterations  required  at  this  plateau  versus  the 
diameter  of  the  sphere  in  wavelengths.  This 
number  of  iterations  for  full  characterisation  seems 
to  rise  roughly  linearly  with  problem  size.  As  the 
matrix  size  varies  with  the  square  of  the  body  size 
or  frequency,  the  cost  of  each  iteration  scales  with 
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the  fourth  power  of  frequency,  giving  a  total  cost 
scaling  in  this  particular  case  with  about  the  fifth 
power  of  frequency. 

4.  Discussion  and  Conclusions 

A  very  simple  modification  of  the  GCR  approach 
has  been  shown  to  be  effective  for  analysis  of 
multiple  right  hand  sides  for  the  large,  dense  and 
unsymmetric  matrices  of  multi-wavelength 
monostatic  res  calculations. 

It  seems  possible  to  obtain  solutions  for  a  large 
number  of  different  illumination  angles  for  a  modest 
multiple  of  the  cost  for  a  single  right  hand  side. 
Typically,  for  the  cases  examined,  this  multiple  is 
~<10  for  an  essentially  unlimited  number  of 
illumination  angles. 

As  noted,  the  cost  scaling  of  iterative  solutions  for 
single  look  angles  is  unclear.  Starting  from  this 
point,  we  conclude  from  the  present  study: 

-  Costs  for  multiple  look  angles  are  independent  of 
the  number  of  look  angles  once  more  than  some 
threshold  number  of  look  angles  is  considered. 

-  This  threshold  seems  to  be  related  to  the  number 
required  to  characterise  the  surface  field,  not  the 
number  required  to  characterise  the  res. 

-  How  the  number  required  to  characterise  the 
surface  field  varies  with  frequency  is  obviously 
geometry-dependant. 

-  On  many  geometries  (eg  the  almond  studied)  the 
surface  field  needs  far  fewer  angles  to  characterise 
it  than  does  the  res. 

-  The  net  result  is  that,  for  a  body  where  this  is 
true,  monostatic  res  characterisation  can  be  obtained 
for  a  small  fraction  of  the  cost  of  repeated  ab  initio 
iterative  solution  of  the  matrix  equation. 
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Figure  1 

NASA  almond,  7  GHz  (6  wavelengths  long);  Iterations  required  versus  number  of 
illumination  angles,  with  angular  range  spanned  by  the  illumination  angles,  measured  from 
head-on,  as  a  parameter. 


Figure  2 
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NASA  almond  (7  GHz,  6  wavelengths  long):  Real  part  of  y -component  of  surface  H  field  at 
surface  location  (-4.1",  0,  0)  as  a  function  of  monostatic  incident  illumination  angle.  The 
illumination  ranges  from  0  to  180°  in  the  equatorial  (VV)  plane  (solid  line),  with  results 
every  6  degrees  (30  uniformly  spaced  illumination  angles)  additionally  marked  by  circles. 
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Figure  3 

NASA  almond  (7  GHz,  6  wavelengths  long):  Monostatic  res  as  a  function  of  monostatic 
incident  illumination  angle.  The  illumination  ranges  from  0  to  180p  in  the  equatorial  (VV) 
plane  (solid  line),  with  results  every  6  degrees  (30  uniformly  spaced  illumination  angles) 
additionally  marked  by  solid  circles.  (Also  shown  are  measured  results,  marked  by  open 
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15  wavelength  dipole:  Iterations  required  versus  number  of  illumination  angles  (distributed 
over  ISCf  in  all  cases).  Inset;  Monostatic  res. 
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Figure  5 

NASA  almond  (7  GHz,  6  wavelengths  long):  Iterations  required  versus  number  of 
illumination  angles,  with  uniformly  spaced  illumination  angles  over  the  range  0  -  16° 
horizontally  and  vertically  from  head-on.  Inset:  Monostatic  (W)  res. 
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Abstract 

The  surface-impedance  boundary  condition  for  the 
Finite-Difference ,  Tune-Domain  (FDTD)  method  is  re¬ 
formulated  using  digital  filtering  theory  and  Z  trans¬ 
forms.  The  approach  expands  upon  recent  work  in  de¬ 
veloping  an  efficient  surface-impedance  boundary  condi¬ 
tion  for  FDTD.  The  present  work  involves  formulating  the 
surface-impedance  boundary  condition  in  the  frequency 
domain  for  a  lossy  dielectric  half-space  and  for  a  thin 
lossy  dielectric  layer  backed  by  a  perfect  conductor.  The 
impedance  function  of  the  lossy  medium  is  approximated 
with  a  series  of  low-pass  filters.  This  approximation  is 
independent  of  material  properties  and  these  low-pass  fil¬ 
ters  are  converted  to  corresponding  digital  filters  using  Z- 
transform  theory.  The  FDTD  surface-impedance  bound¬ 
ary  condition  is  reformulated  in  the  Z  domain ,  and  the 
corresponding  time-domain  electric  field  sequence  updat¬ 
ing  equation  involves  a  recursive  formula.  Results  are 
presented  for  both  one  and  two-dimensional  test  prob¬ 
lems. 


1  Introduction 

The  surface-impedance  boundary  condition  (SIBC)  is 
a  frequency-domain  concept  that  is  used  to  simplify  scat¬ 
tering  calculations  by  eliminating  the  internal  volume  of 
lossy  dielectric  objects.  The  SIBC  relates  tangential  elec¬ 
tric  and  magnetic  field  components  at  the  surface  of  the 
object  through  an  impedance  which  is  a  function  of  the 
object* s  material  parameters.  Thus,  the  material  proper¬ 
ties  of  the  object  are  accounted  for,  and  if  exterior  fields 
are  only  of  interest,  the  SIBC  eliminates  the  need  to  cal¬ 
culate  fields  inside  the  scatterer.  This  results  in  a  savings 
in  computational  resources  by  reducing  the  computational 
storage  requirements  and/or  computation  time. 

To  analyze  electromagnetic  field  interaction  with  lossy 
dielectric  objects,  the  Finite-Difference,  Time-Domain 
(FDTD)  method  requires  that  the  interior  of  the  object  be 


modeled  in  order  for  fields  to  penetrate  the  body.  Since 
the  wavelength  inside  a  lossy  dielectric  material  is  much 
less  than  the  free  space  wavelength,  accurate  modeling 
often  requires  a  fine  spatial  grid  resulting  in  a  relatively 
large  number  of  cells  for  moderately  sized  objects.  For 
calculations  where  only  exterior  response  is  of  interest,  a 
conducting  dielectric  object  can  be  replaced  by  a  SIBC 
over  the  surface  of  the  object.  Thus,  this  boundary  condi¬ 
tion  eliminates  the  spatial  quantization  of  the  object  and 
reduces  the  overall  size  of  the  solution  space  by  elimi¬ 
nating  cells  within  the  object  and  by  allowing  larger  cells 
to  be  used  in  the  exterior  region.  The  larger  cells  reduce 
the  storage  requirements  since  fewer  cells  are  required  to 
model  the  same  physical  dimensions  of  the  object.  The 
computation  time  for  the  reduced  solution  space  is  also 
decreased  because  fields  in  cells  within  the  object  are  not 
updated. 

Most  applications  of  the  SIBC  have  traditionally  been 
in  the  frequency  domain  [1]— [18].  Recently,  time-domain 
surface-impedance  concepts  have  received  considerable 
attention  in  the  literature  [19].  There  have  been  several 
FDTD  implementations  of  the  surface-impedance  bound¬ 
ary  condition  introduced  by  various  authors  [20]— [38]. 
Each  implementation  has.  certain  advantages  and  disad¬ 
vantages,  but  all  strive  to  obtain  the  most  efficient  and  ac¬ 
curate  method.  These  FDTD  surface-impedance  bound¬ 
ary  conditions  are  reviewed  in  a  separate  article  [39]. 

It  has  been  demonstrated  recently  that  dispersive  and 
nonlinear  optical  media  can  be  modeled  in  FDTD  using 
digital  filtering  and  Z-transform  theory  [40].  Materials 
with  Debye  or  Lorentz  dispersion  have  rational  functions 
of  frequency  for  the  dielectric  permittivity.  These  func¬ 
tions  that  define  relative  permittivities  as  a  function  of 
frequency  have  direct  Z-transforms,  thus  allowing  the  re¬ 
lationship  between  electric  flux  density  and  field  intensity 
to  be  formulated  directly  in  the  Z-domain.  The  result¬ 
ing  time-domain  updating  equations  for  the  sequences  in¬ 
volve  recursive  formulas  and  are  computationally  identi¬ 
cal  to  the  time-domain,  differential-equation  based  meth- 
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ods.  Excellent  results  have  been  obtained  in  modeling  De¬ 
bye,  Lorentz  and  nonlinear  optical  media.  The  difficulty 
in  applying  Z-transform  theory  to  the  surface-impedance 
boundary  condition  is  that  the  frequency-domain,  surface- 
impedance  function  is  irrational.  With  the  recent  work 
of  Oh  et  al  [37],  the  normalized  irrational  surface- 
impedance  function  is  approximated  with  a  series  of  first- 
order,  rational  functions.  This  approximation  is  inde¬ 
pendent  of  material  properties,  and  the  first-order,  ratio¬ 
nal  functions  can  be  transformed  into  an  equivalent  Z- 
domain  form.  This  paper  extends  the  work  in  [37]  by 
providing  the  technical  approach  for  reformulating  the 
surface-impedance  boundary  condition  using  digital  fil¬ 
tering  and  Z-transform  theory.  The  surface-impedance 
boundary  condition  is  converted  to  the  Z-domain  by  two 
methods  and  the  resulting  recursive  updating  formulas  for 
the  electric  field  intensity  sequence  are  almost  computa¬ 
tionally  identical  to  that  in  [37].  There  are  several  advan¬ 
tages  of  using  the  Z-transform  approach  to  implement  the 
SIBC  for  the  FDTD  method.  The  Z-transform  approach 
provides  a  more  accurate,  intuitive  and  consistent  method 
to  implement  the  SIBC  based  upon  the  discrete  nature  of 
the  FDTD  method.  In  the  traditional  recursive  convolu¬ 
tion  approach,  the  SIBC  is  implemented  by  convolving 
analog  time-domain  impedance  and  magnetic  field  signals 
that  have  been  discretized  in  time.  However,  the  discrete 
nature  of  the  FDTD  method  results  in  discrete  impedance 
and  magnetic  field  sequences  which  are  manipulated  eas¬ 
ily  and  accurately  using  Z-transforms.  The  Z-domain 
SIBC  system  functions  are  digital  filters  that  operate  on 
the  discrete  FDTD  magnetic  field  sequences  which  pro¬ 
vides  a  more  accurate  implementation  and  a  more  intu¬ 
itive  and  cohesive  theoretical  formulation. 

In  this  paper.  Section  2  develops  the  surface  impedance 
boundary  condition  in  the  Z-domain  and  provides  the  re¬ 
cursive  updating  equations  for  the  electric  field  inten¬ 
sity  sequence  for  both  a  lossy  dielectric  half- space  and 
a  thin  lossy  dielectric  layer  backed  by  a  perfect  electri¬ 
cal  conductor  (PEC).  Section  3  presents  one-  and  two- 
dimensional  results,  and  Section  4  provides  concluding 
remarks. 


2  Z  Transform  Approach 

This  section  develops  the  Z-transform  approach  for 
both  a  lossy  dielectric  half-space  and  a  thin,  PEC-backed, 
lossy-dielectric  layer.  The  planar  first-order,  frequency- 
domain  SIBC  is  used  and  then  an  efficient  implementation 
of  the  corresponding  time-domain  SIBC  is  developed  us¬ 
ing  digital  filtering  and  Z-transform  theory.  This  approach 
is  an  extension  of  previous  work  by  Oh  et  al.  [37]  in  de¬ 
veloping  an  efficient  time-domain  SIBC  using  recursion. 
The  notation  presented  in  that  paper  [37]  has  been  pre¬ 


served  as  much  as  possible  in  the  present  work. 


2.1  Lossy  Dielectric  Half-Space 

The  first-order  (or  Leontovich)  impedance  boundary 
condition  relates  tangential  total  field  components  and  is 
given  in  the  frequency  domain  as  [4] 

E(oj)  -  n  [n  •  J5(w)|  =  Zt(u)  [n  x  (1) 

where  u>  is  the  radian  frequency,  n  is  the  unit  outward  nor¬ 
mal  from  the  surface  and  Zs{u)  is  the  frequency-domain 
surface-impedance  of  the  material.  An  ejut  time  depen¬ 
dence  is  assumed  and  suppressed.  This  frequency-domain 
SIBC  is  for  a  planar  material  interface  and  does  not  ac¬ 
count  for  the  surface  curvature  of  an  object.  Since  the 
present  formulation  uses  the  planar  SIBC,  it  is  limited 
to  those  geometries  where  the  smallest  radius  of  curva¬ 
ture  is  relatively  large  compared  to  the  wavelength.  The 
frequency-domain  surface-impedance  is  given  by 


Zs(u)  — 


I  juy 
a  +ju>e 


(2) 


When  the  impedance  boundary  is  between  free  space  and 
the  dielectric,  it  is  assumed  that  the  complex  refractive 
index,  N,  obeys  the  restriction 

|iV|  »  1  (3) 


where  N  =  vVrer(w).  This  restriction  is  imposed  so 
that  the  wave  impedance  in  the  material  is  independent 
of  the  incidence  angle  and  will  be  equal  for  both  polar¬ 
izations.  This  restriction  limits  the  applicability  of  the 
present  formulation  to  those  media  with  large  conductiv¬ 
ity  or  permittivity,  but  this  restriction  is  feasible  for  most 
practical  simulations  where  a  SIBC  would  be  applicable. 
For  cases  of  low  \N\,  the  approach  developed  by  Kellali 
et  al.  [29]  includes  the  incidence  angle  in  the  SIBC.  That 
formulation  may  be  beneficial  in  cases  of  lower  |iV|  to 
avoid  gridding  and  updating  fields  within  a  large  object 
Since  the  restriction  of  (3)  holds,  then  (2)  is  applicable 
for  both  polarizations  in  the  two-  and  three-dimensional 
cases.  This  implies  that  the  transmitted  fields  inside  the 
object  will  propagate  almost  normal  to  the  object  surface; 
hence  the  wave  impedance  will  be  equal  for  both  polar¬ 
izations. 

The  frequency-domain  surface-impedance  function  of 
(2)  can  be  rewritten  using  the  complex  frequency  variable 
s  =  ju>  as 

where  a  =  a/e  and  Zt  is  the  intrinsic  wave  impedance  of 
the  material  given  by  Z.  =  \j \ij  Following  the  proce¬ 
dure  in  [37],  equation  (4)  can  be  rewritten  as  a  normalized 
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impedance  function  given  by 


Zn(s') 


Zs(s') 

Zi 


(5) 


where  s'  =  s/a.  The  normalized  surface-impedance 
function  is  approximated  in  the  frequency  domain  by  a 
series  of  first  order  rational  functions  of  the  form 


L 

ZnW)*  l-£ 


1  =  1 


Cl 

s'  -f  u ;/ 


(6) 


This  approximation  is  over  the  real  axis  interval  sl  ~ 
[0,3],  which  will  accommodate  most  materials  up  to  sev¬ 
eral  tens  of  Gigahertz.  The  residues  Ci  and  poles  cj/ 
are  given  for  sixth,  seventh  and  eighth-order  approxi¬ 
mations  in  [37].  Figure  1  shows  the  percentage  error 
for  the  eighth-order  approximation,  which  will  be  used 
throughout  this  paper.  Using  the  approximation  in  (6), 


representation  of  each  first-order  rational  function  is  a  de¬ 
caying  exponential.  As  a  result,  the  convolution  with  the 
tangential  magnetic  field  can  be  updated  using  recursion 
to  avoid  storing  the  complete  time  history  of  the  magnetic 
field  at  the  material  interface.  This  recursive  approach 
was  first  proposed  in  [22]-[23]  and  has  been  expanded 
upon  by  several  authors  [29]— [31],  [34]— [38].  Sullivan 
[41]  has  shown  the  multiplication  theorem  for  convolu¬ 
tion  involving  FDTD  discrete  field  sequences  includes  the 
sampling  factor  T.  Thus  the  Z-domain  surface-impedance 
boundary  condition  is  given  by 


Et{z)  =  Zs(z)  [nxff(z)]T 


(9) 


Each  of  the  first-order  rational  functions  in  (7)  is  a  low- 
pass  filter.  When  an  analog  filter  function  is  known,  dig¬ 
ital  filters  can  be  obtained  directly  using  three  different 
approaches  [42].  These  approaches  are  discussed  in  the 
following  sections. 


Figure  1:  Percent  error  versus  complex  frequency  (s')  in 
eighth-order  approximation  to  surface  impedance  func¬ 
tion. 


the  frequency-domain,  surface-impedance  function  can 
be  written  as 


8 


1=1 


aZiCi 
s  +  aui 


(7) 


2.1.1  Impulse  Invariance  Method 

The  first  approach  to  obtaining  a  digital  filter  from 
an  analog  filter  is  to  use  the  impulse-invariance  proce¬ 
dure  which  involves  choosing  the  unit  sample  response 
of  the  digital  filter  as  equally  spaced  samples  of  the  im¬ 
pulse  response  of  the  analog  filter.  The  time-domain, 
surface-impedance  impulse  response  is  found  via  an  in¬ 
verse  Laplace  transform  of  the  frequency-domain  expres¬ 
sion  in  (7)  to  be 


8 

zs(t)  =  Zi5(t)  -  J2  aCie-auit  (10) 

1=1 

where  6(t)  is  the  Dirac  delta  function.  If  this  time-domain 
impulse  response  is  now  sampled  with  interval  T,  the  dis¬ 
crete  surface-impedance  impulse  response  is 

8 

zs(n)  =  z$(nT)  =  Zi5{nT)  -  ]T  aCie~au'nT  (11) 

1=1 

Taking  the  Z-transform  of  both  sides  of  this  equation  gives 
the  Z-domain,  surface-impedance  function  of 


The  time-domain,  surface-impedance  boundary  condition 
is  obtained  by  applying  the  Convolution  Theorem  to  (1) 
which  results  in 


ry  (  \  _ 

s[Z)  -  T  2_.  !  _  Z-1  e-autT 


(12) 


et{t)  =  zs(t)  ®  [n  x  £(t)j  (8) 

where  the  0  is  the  convolution  operator,  e(t)  and  h(t)  are 
the  time-domain  electric  and  magnetic  field  intensities  at 
the  impedance  boundary  and  the  subscript  denotes  tan¬ 
gential  field.  From  (7),  it  is  clear  that  the  time-domain 


The  Z-domain  surface-impedance  boundary  condition  of 
(9)  is  then 


Et(z)  = 


*-Er 


aQT 


1=1 


z-le-auJiT 


n  x  JT(z)j 


(13) 
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This  equation  can  be  rewritten  after  some  algebra  as 

8 

Et  (z)  =  Zi\nxH  (z)  1  (2)  (14) 

L  Z= 1 

where 

Fi(z)  =  e~au,T z~1  Fi(z)  +  aQZi  [n  x  H(z)]  (15) 

is  the  recursive  updating  equation  in  the  Z-domain.  Re¬ 
call  that  the  z-1  term  is  a  delay  operator,  and  in  the  time 
domain,  these  equations  become 

e?(l fc)  =  Zi  \n  x  hn~l/2(k  -  1/2)1  -  Z (16) 

1  1=1 

with 

fP(k)  =  e~aLJlTfi~\k)  (17) 

+aCtZi  [n  x  hn~1/2(k  -  1/2)] 

as  the  recursive  updating  equation  suitable  for  FDTD  im¬ 
plementation.  Note  that  equation  (16)  is  identical  to  equa¬ 
tion  (8a)  in  [37]  with  a  minor  change  in  notation.  By  us¬ 
ing  the  impulse  invariance  procedure,  it  is  assumed  that 
the  field  sequences  are  piecewise  constant  in  time.  This 
will  result  in  only  a  first-order  accurate  boundary  con¬ 
dition,  similar  to  the  recursive  convolution  approach  for 
dispersive  media  involving  piecewise  constant  field  com¬ 
ponents.  When  implemented  and  tested,  the  impulse  in¬ 
variance  procedure  led  to  instabilities  because  of  the  alias¬ 
ing  problems  in  the  digital  filter  frequency  response,  and 
because  the  time  step  chosen  for  the  FDTD  calculations 
was  not  small  enough  to  resolve  the  characteristic  time 
constants  of  some  of  the  exponential  terms  in  the  unit  im¬ 
pulse  response.  It  is  interesting  to  note  that  the  impulse- 
invariance  design  approach  was  shown  to  be  unstable  for 
FDTD  modeling  of  dispersive  media  for  certain  materi¬ 
als  exhibiting  Lorentz  dispersion  [43].  Therefore,  no  re¬ 
sults  are  shown  using  the  impulse-invariance  procedure, 
but  this  section  was  included  for  completeness  and  for 
comparison  to  the  other  Z-domain  methods. 


2.1.2  Backward  Difference  Method 

With  the  frequency-domain  surface-impedance  of  (7), 
the  surface-impedance  boundary  condition  can  be  rewrit¬ 
ten  to  be 


Et(s) 


Zi-E 


i=i 


aCiZi 
s  +  a  ui 


(18) 


The  backward  difference  method  is  a  digital  filter  design 
technique  based  upon  the  numerical  solution  of  a  differen¬ 
tial  equation.  Each  of  the  first-order  rational  terms  above 


can  be  thought  of  as  the  analog  system  function  for  a  first- 
order  time-domain  differential  equation.  By  approximat¬ 
ing  the  time  derivatives  in  the  differential  equation  by  a 
backward  difference,  the  digital  system  function  is  ob¬ 
tained  from  the  analog  system  function  by  a  substitution 
of  variables  . 


Substituting  this  into  (18)  and  after  some  manipulation 
gives  the  Z-domain  surface-impedance  boundary  condi¬ 
tion 


Et(z) 


aQZiT 

Zi~f^  (l-z-^  +  aunT 
[n  x  tf(z)] 


(20) 


Note  that  an  extra  T  term  does  not  appear  above  because 
this  transformation  does  not  involve  the  convolution  theo¬ 
rem,  hence  the  T  factor  is  not  present.  After  some  algebra, 
equation  (20)  can  be  written  in  the  form 


Et  (z)  =  Zi  In  xH (z)]  (z)  (21) 

L  Z= 1 


where 


Fi(z)  = 


1  +  awiT 


Ft(z)  + 


aCiZjT 

1  +  auiT 


[nxH(z)}  (22) 


is  the  recursive  updating  equation  in  the  Z-domain.  The 
time-domain  SIBC  is  the  same  as  that  given  in  (16)  with 


f?(k)  = 


as  the  recursive 
plementation. 


1 


/r  W 


(23) 


1  4-  au>iT 
aCiZ{T 
+ 1  +  aioiT 

updating  equation  suitable  for  FDTD  im- 


fn  x  hn~1/2{k  -  1/2)] 


2.1.3  Bilinear  Transformation  Method 


If  the  time-domain  differential  equation  used  for  the 
backward  difference  method  for  each  of  the  first-order  ra¬ 
tional  functions  in  (18)  is  integrated,  and  then  this  integral 
is  approximated  by  the  trapezoidal  rule,  the  corresponding 
Z-domain  system  function  can  be  obtained  from  the  ana¬ 
log  function  by  a  different  substitution  of  variables.  This 
is  the  bilinear  transformation  method,  and  the  substitution 
of  variables  is  . 


Performing  this  substitution  of  variables  in  (18)  and  af¬ 
ter  some  algebra  gives  the  Z-domain  surface-impedance 
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boundary  condition  to  be 


Et(z)  «  [Zt- 
®  aQZiT  (1  +  z-1) 

^  (2  +  awiT)  -  (2  -  auiT)  z"1 


(25) 

x  i?(z)| 


Note  that  an  extra  T  term  does  not  appear  above  because 
the  bilinear  transformation  does  not  involve  the  convolu¬ 
tion  theorem,  hence  the  T  factor  is  not  required.  Equation 
(25)  can  be  written  in  the  same  form  as  (14)  where 


is  the  recursive  updating  equation  in  the  Z-domain.  The 
corresponding  time-domain  SIBC  update  equation  is  the 
same  as  (16)  and  the  recursive  time-domain  update  equa¬ 
tion  is  given  by 


+(ra) 

+hn~3/2(k  -  1/2))]  (27) 

and  is  suitable  for  FDTD  implementation.  In  the  paper 
by  Oh  et  al ,  the  recursive  updating  equation  is  given  by 
equation  (8b)  of  that  paper  [37]  (with  a  slight  change  of 
notation)  as 


8 

«?(*)  =  Zi  [n  x  hn-l>2{k  -  1/2)]  -  Y,  Ai(k)  (28) 


and  stored  to  maximize  efficiency.  Comparing  the  coeffi¬ 
cients  in  (17),  (23)  with  (29),  it  is  clear  that  the  former  two 
equations  are  first-order  accurate  since  they  only  involve 
one  past  time  level  of  the  tangential  magnetic  field.  Note 
that  the  bilinear  transformation  recursion  equation  in  (27) 
requires  one  less  multiply  and  one  less  coefficient  storage 
location  per  recursive  updating  variable  than  the  equation 
in  (29).  This  makes  the  bilinear  transformation  approach 
slightly  more  efficient  than  (29).  The  bilinear  transforma¬ 
tion  maps  the  entire  left  half  of  the  complex  plane  inside 
the  unit  circle  and  the  imaginary  axis  in  the  complex  plane 
becomes  the  unit  circle.  As  a  result,  frequency  warping 
will  take  place  when  transferring  from  the  analog  system 
to  the  digital  system.  The  bilinear  transformation  is  most 
useful  when  this  distortion  can  be  tolerated  or  compen¬ 
sated.  When  designing  a  digital  filter  using  the  bilinear 
transformation,  the  analog  cutoff  frequencies  can  be  pre¬ 
warped  so  that  the  digital  cutoff  frequencies  will  fall  at 
the  correct  design  cutoff  frequencies.  However,  in  this 
application,  this  distortion  in  the  frequency  axis  can  be 
tolerated  and  no  compensation  is  required. 

2*2  PEC-Backed  Thin  Lossy  Dielectric 
Layer 

Following  the  notation  in  [37],  this  section  devel¬ 
ops  a  Z-transform  SIBC  implementation  for  a  PEC- 
backed,  thin,  lossy  dielectric  layer.  The  geometry  is  one¬ 
dimensional  and  the  lossy  dielectric  layer  has  thickness, 
d,  and  parameters  e*,  fxi  and  at  #  0.  It  is  assumed  that  the 
intrinsic  impedance  of  the  dielectric  layer,  Zi(u)7  (equa¬ 
tion  (2))  obeys  the  restriction  given  in  (3)  and  that  d  is 
less  than  one-half  the  cell  size.  Since  the  dielectric  layer 
is  backed  by  a  PEC,  the  surface  impedance  looking  into 
the  layer  is 


where  Af(k)  is  the  recursive  updating  variable  given  by 
Af(k)  =  ^  [1  +  (e~aUiT  -  1)  /(auiT)]  ■ 

UJi 

[n  x  hn~V'\k-  1/2)]  + 

—  [1  ftauiT)  -  e~aUiT(l  +  1  /(auiT))]  ■ 

u>i 

\n  x  hn~3'2{k  -  1  /2)]  +  e~aUiTA^~l  ( k )  (29) 

The  coefficients  of  this  recursive  updating  equation  were 
obtained  by  directly  evaluating  the  convolution  integral 
of  the  time-domain  surface-impedance  function  with  the 
tangential  magnetic  field  next  to  the  impedance  bound¬ 
ary.  The  coefficients  in  (17),  (23)  and  (27)  were  obtained 
by  manipulating  the  surface-impedance  function  in  the  Z- 
domain.  All  of  the  coefficients  in  the  recursive  updating 
equations  of  (17),  (23),  (27)  and  (29)  can  be  precomputed 


Zs(u)  =  Zi(u)  tan(7(cj)d)  (30) 

where  y(u)  is  the  propagation  constant  If  d  is  small,  then 
the  approximation 

tan(7(cj)d)  ^7 (u)d  (31) 

is  applied  to  (30)  to  give 

Z$(cu )  =  Zi(cj)j(cj)d  =  jufxid  (32) 

Using  the  Laplace  transform  variable  s  =  ju  gives 

Zs(s)  —  fiids  (33) 

Now  that  the  thin-layer  surface-impedance  is  expressed 
in  terms  of  s,  the  Z-transform  methods  can  be  applied  di¬ 
rectly  to  obtain  the  corresponding  update  equations. 
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2.2.1  Backward  Difference  Method 

For  the  backward  difference  method,  the  substitution  of 
variables  given  in  (19)  is  used  in  (33)  to  give 

ZM=„i(±=f^)  (34) 

in  the  Z-domain.  Now  the  Z-domain  SIBC  can  be  written 
as 

Et(z)  =  ^  [(1  -  z~l)  (n*H (z))]  (35) 

and  the  time-domain  SIBC  is  given  by 

e7  (Jfe)  =  ^[nx(H2M/2) 

-hn~3/*(k  -  1/2))]  (36) 

which  is  suitable  for  FDTD  implementation.  Note  that 
this  update  equation  is  identical  to  (12)  in  [37].  The  back¬ 
ward  difference  method  is  the  only  stable  Z-transform 
method  for  the  PEC-backed,  lossy  dielectric  layer.  The  bi¬ 
linear  transformation  method  was  unstable  because  the  re¬ 
sulting  Z-domain  surface-impedance  function  has  a  pole 
at  =  7 r  where  fl  is  the  digital  filter  frequency  variable. 

3  Demonstration 

3.1  Lossy  Dielectric  Half-Space 


frequency  for  a  pulse  normally  incident  on  a  lossy  dielec¬ 
tric  half  space.  The  problem  space  size  was  100  cells, 
the  impedance  boundary  was  located  at  cell  100  and  the 
electric  field  was  sampled  at  cell  100.  The  maximum  fre¬ 
quency  of  interest  was  10  GHz  and  the  incident  electric 
field  was  a  Gaussian  pulse  with  unity  amplitude  of  the 
form  . 

Ei(t)  =  e-^~to)/T)  u(t)  (37) 

with  to  =  80  ps  and  r  =  20  ps.  The  pulse  was  win¬ 
dowed  in  time  at  approximately  -80  dB  with  a  rectan¬ 
gular  window  width  of  64  time  steps.  The  frequency 
response  of  the  pulse  contained  significant  energy  to  12 
GHz.  Two  computations  were  made  with  a  =  2.0  S/m  and 
a  =  200.0  S/m  corresponding  to  loss  tangents  at  10  GHz 
of  3.6  and  360,  respectively.  The  permittivity  and  per¬ 
meability  of  the  lossy  dielectric  half  space  were  taken  as 
those  of  free  space.  The  cell  size  and  time  step  were  750 
/xm  (40  cells/ Ao  at  10  GHz)  and  2.5  ps  respectively,  and  a 
total  of  8192  time  steps  were  evaluated.  For  each  FDTD 
computation,  a  reflection  coefficient  versus  frequency  was 
obtained  by  dividing  the  Fourier  transform  of  the  scat¬ 
tered  field  by  the  Fourier  transform  of  the  incident  field 
at  cell  100.  The  incident  field  was  obtained  by  running 
the  FDTD  code  with  free  space  only  and  recording  the 
electric  field  at  cell  100.  The  scattered  field  is  then  ob¬ 
tained  by  subtracting  the  time-domain  incident  field  from 
the  time-domain  total  field.  The  results  are  compared  with 
the  analytic  surface-impedance  reflection  coefficient  com¬ 
puted  from 


To  demonstrate  these  different  approaches,  equation 
(16),  along  with  the  recursive  updating  equations,  (23) 
and  (27)  were  implemented  in  a  one-dimensional  total 
field  FDTD  code  for  the  geometry  shown  in  Figure  2.  The 
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Figure  2:  Problem  geometry  for  FDTD  SIBC  with 
one-dimensional  FDTD  grid  and  the  placement  of  the 
impedance  boundary. 

problem  was  to  calculate  the  reflection  coefficient  versus 


\R\  = 


1  -  -y/1  -  jtr/^e o 
1  +  y/l  -  jcr/weo 


(38) 


where  a  is  the  conductivity  of  the  dielectric  half-space. 
The  Z-transform  results  are  also  compared  with  the  di¬ 
rect  recursive  convolution  approach  of  Oh  et  al.  [37] 
using  the  eighth  order  approximation  and  with  a  conven¬ 
tional  FDTD  calculation  for  a  simulated  lossy  half-space. 
In  all  figures,  the  following  abbreviations  apply:  recur¬ 
sive  convolution,  RC;  backward  difference,  BD;  bilinear 
transformation,  BLT.  Figures  3  and  4  show  the  SIBC  re¬ 
flection  coefficient  results  for  all  methods  compared  with 
the  analytic  SIBC  result  for  a  —  2.0  S/m  and  cr  =  200.0 
S/m,  respectively.  Notice  the  agreement  is  excellent  for 
the  o-  =  2.0  case  and  is  less  at  higher  frequencies  in  the  o' 
=  200.0  case.  This  discrepancy  is  due  to  discretization  er¬ 
ror.  In  Figure  4,  note  that  the  SIBC  implementations  have 
approximately  the  same  absolute  accuracy  as  the  FDTD 
result  and  the  SIBC  cannot  be  expected  to  perform  any 
better  than  the  FDTD  calculation.  Therefore,  these  re¬ 
sults  are  encouraging.  Comparing  Figure  4  to  Figure  3c 
in  [37],  it  is  clear  that  Figure  3c  exhibits  better  agree¬ 
ment  at  higher  frequencies.  However,  this  is  because  the 
reflection  coefficient  shown  in  that  graph  was  computed 
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Figure  3:  Reflection  coefficient  magnitude  versus  fre¬ 
quency  for  normal  incidence  plane  wave  reflection  from 
a  lossy  dielectric  half  space  with  a  =  2.0  S/m  using  the 
analytic  andFDTD  SIBC  methods. 


from  a  closed  form  representation  for  the  rational  func¬ 
tion  approximation.  Note  in  Figure  4  from  [37]  that  the 
results  for  the  case  of  a  —  2.0  S/m  with  a  six-term  ratio¬ 
nal  function  approximation  exhibit  the  same  type  of  high 
frequency  behavior  as  the  results  presented  in  Figure  3 
and  Figure  4.  Thus,  the  results  presented  here  are  consis¬ 
tent  with  those  presented  in  [37].  The  Z- transform  SIBCs 
are  only  first-order  accurate  overall  because  the  magnetic 
field  at  the  impedance  boundary  is  approximated  by  the 
magnetic  field  one-half  cell  in  front  of  the  impedance 
boundary.  Although  it  was  anticipated  that  the  bilinear 
transformation  method  was  second-order  accurate  in  time, 
when  applied  in  practice,  it  is  only  first-order  accurate. 
This  is  illustrated  in  Figure  5  where  the  BLT  method  is 
used  at  a  cell  spacing  of  325  pm  (1/2  the  previous  cell 
size)  with  a  conductivity  of  a  -  200.0  S/m.  Note  that 
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Figure  4:  Reflection  coefficient  magnitude  versus  fre¬ 
quency  for  normal  incidence  plane  wave  reflection  from 
a  lossy  dielectric  half  space  with  a  -  200.0  S/m  using  the 
analytic  and  FDTD  SIBC  methods. 


Figure  5:  Reflection  coefficient  magnitude  versus  fre¬ 
quency  for  normal  incidence  plane  wave  reflection  from 
a  lossy  dielectric  half  space  with  a  =  200.0  S/m  using  the 
FDTD  BLT  SIBC  at  grid  resolutions  of  Sz  and  5zf2. 

the  new  result  is  approximately  50%  closer  to  the  analyt¬ 
ical  solution.  Similar  behavior  was  observed  with  the  BD 
method. 

3.2  PEC-Backed  Thin  Lossy  Dielectric 
Layer 

The  calculation  parameters  for  the  thin  layer  example 
are  the  same  as  for  the  previous  example  except  the  mate¬ 
rial  parameters  have  the  following  values: 


ei  =  2eo 

Pi  —  2po 

u\  —  2.0  S/m  (39) 


and  d  =  0.15z.  In  this  example,  the  impedance  is  calcu¬ 
lated  at  the  thin-layer  interface  by  calculating  the  ratio  of 
the  electric  and  magnetic  fields  computed  via  the  SIBC. 
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This  is  then  compared  with  the  exact  impedance  value 
computed  from  (30).  Excellent  agreement  is  obtained  for 
the  BD  Z-transform  method  as  shown  in  Figure  6. 
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Figure  6:  Magnitude  of  complex  impedance  for  a  thin, 
PEC-backed  lossy  dielectric  layer  with  e  =  2e0,  —  2/jq, 

a  =  2.0  S/m  and  d  =  0.1 5z. 


3.3  Two-Dimensional  Scattering  Example 

In  this  example,  a  two-dimensional  TM  FDTD  code 
was  adapted  to  use  the  recursive  convolution  and  Z- 
transfoim  SIBCs.  Figure  7  shows  the  two-dimensional 
geometry  for  both  the  plane  wave  and  point  source  ex¬ 
citations.  A  plane  wave  is  incident  from  the  —  x  axis 

I  y 


the  following  parameters:  er  =  1,  Hr  =  1  and  a  =  20.0 
S/m.  The  problem  space  size  was  300  by  300  cells,  the 
cell  size  is  500  /mi  and  the  time  step  is  1.2  ps.  Scattering 
angles  of  <p  =  180,  90  and  0  degrees  were  used,  where  <j> 
is  measured  from  the  +x  axis.  The  incident  pulse  was  of 
the  same  form  as  (37)  with  to  =  37.7  ps  and  r  =  9.45 
ps  and  a  second-order  Liao  absorbing  boundary  condition 
was  used  [44].  Since  an  analytical  solution  is  not  avail¬ 
able,  each  Z-transform  approach  is  compared  with  a  con¬ 
ventional  FDTD  calculation  for  the  same  geometry.  Fig¬ 
ure  8  shows  the  monostatic  scattering  width  at  <p  =  180° 
and  Figures  9  and  10  show  the  bistatic  scattering  width 
at  <t>  =  90°  and  <f>  =  0°,  respectively.  Note  the  agree- 
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Figure  8:  Monostatic  scattering  width  at  <j>  =  180°  for 
lossy  square  cylinder  using  conventional  FDTD  and  the 
various  SIBC  implementations. 
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Figure  7:  Problem  geometry  for  two-dimensional  TM 
scattering  calculations  showing  square  cylinder,  incident 
plane  wave,  point  source  and  field  sampling  points. 

and  a  two-dimensional  scattering  width  is  calculated  for 
a  square  cylinder.  The  cylinder  is  10  cm  square  and  has 
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Figure  9:  Bistatic  scattering  width  at  <j>  =  90°  for  lossy 
square  cylinder  using  conventional  FDTD  and  the  various 
SIBC  implementations. 

ment  is  good  in  all  three  cases.  It  is  interesting  to  note 


21 


Figure  10:  Bistatic  scattering  width  at  <j>  =  0°  for  lossy  Figure  1 1 :  Near  zone  scattered  electric  field  at  a  point  25 
square  cylinder  using  conventional  FDTD  and  the  various  cells  above  a  10  cm  lossy  square  cylinder  with  a  point 
SIBC  implementations.  source  excitation. 


that  the  half-cell  spacing  error  between  electric  and  mag¬ 
netic  fields  in  the  SIBC  seems  to  manifest  itself  more  in 
the  monostatic  scattering  width  of  Figure  7  than  in  the 
bistatic  scattering  width  patterns.  This  is  most  likely  due 
to  cancellation  of  the  SIBC  error  sources  from  all  four 
sides  of  the  square  cylinder  in  the  bistatic  directions.  A 
nonplanar  wave  example  was  also  used  with  the  excita¬ 
tion  point  located  at  grid  point  io,  jo  (see  Fig.  7)  and  a 
soft  excitation  source  of  the  form 

Etotal(i0 ,  jo)  -  E{dtd{i0 ,  jo)  +  Etnc (i0,  jo)  (40) 

where  E(dtd(i0Jo)  is  the  free-space  electric  field  up¬ 
dated  using  FDTD  and  Elnc(io,jo)  is  the  source  term 
of  the  same  form  as  (37).  The  source  point  was  located 
at  io  =  25,  jo  =  151  and  the  pulse  parameters  were 
to  =  37.7  ps  and  r  =  9.45  ps  with  the  amplitude  of 
the  source  at  1000  V/m.  The  results  are  shown  in  Fig¬ 
ure  1 1  for  a  near-field  sampling  point  located  at  grid  point 
(151,276).  Note  the  agreement  is  excellent  for  the  recur¬ 
sive  convolution  method  and  for  both  Z-transform  meth¬ 
ods. 

4  Conclusion 

An  efficient  implementation  of  the  frequency- 
dependent  SIBC  for  the  FDTD  method  based  upon 
Z-transforms  has  been  presented.  Both  the  backward 
difference  and  bilinear  transformation  Z-transform  meth¬ 
ods  were  implemented  and  tested  for  a  lossy  dielectric 
half-space  and  the  backward  difference  method  was 
implemented  and  tested  for  a  thin  PEC -backed  lossy 
dielectric  layer.  Excellent  agreement  was  obtained  on 
one-  and  two-dimensional  TM  problems  using  both 


plane  and  nonplanar  wave  excitations.  Although  only  a 
square  geometry  was  considered  in  the  two-dimensional 
example  in  this  paper,  scattering  from  circular  cylinders 
was  demonstrated  using  a  similar  SIBC  in  [27].  The 
extension  to  the  two-dimensional  TE  case  and  the  three- 
dimensional  case  is  straightforward. 
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Abstract  -  NEC  4  A  has  been  used  to  compute  the  responses  of  the  antennas  in  a  shipboard  high  frequency  direction  finding  system  which 
employs  the  CIDF  algorithm  to  derive  bearing  estimates .  This  paper  discusses  the  computational  results  as  well  as  the  performance  of  the 
simulation  in  which  the  results  were  utilized . 


1 .  INTRODUCTION 
1.1.  Background 

The  requirement  to  carry  out  high 
frequency  direction  finding  (HFDF) 
from  both  shore  stations  and  ships 
arose  during  WWII.  As  recounted  by 
Price  [1]  and  Redgment  [2],  [3],  it 

provided  a  means  for  countering  the 
German  submarine  threat  to  allied 
merchant  ship  convoys.  The  anti¬ 
submarine  mission,  as  well  as  other 
modern  day  missions,  continue  to  make 
shipboard  HFDF  an  important 
consideration  in  the  design  of 
electromagnetic  systems  for  many 
ships. 

Direction  finding  from  a  ship  at  high 
frequencies  (HF)  is  a  challenging 
problem  because  in  the  HF  band, 
antennas  may  interact  strongly  with 
the  ship's  superstructure  and  their 
in-situ  phase  and  amplitude  responses 
can  deviate  significantly  from  their 
free  space  responses.  The  bearing 
errors  observed  using  simple  systems 
can  be  10°-20°  and  in  some  cases, 
error  curves  may  even  be  multivalued 
or  discontinuous,  as  reported  by 
Crampton  et  al.  [4],  and  Travers  et 
al.  [5].  As  a  result,  direction 
finding  (DF)  techniques  which  can  be 
utilized  at  higher  frequencies  do  not 
work  well  aboard  ship  at  HF.  At  HF , 
one  must  use  a  technique  which 
accounts  for  the  interaction  of  the 
antennas  with  the  ship's 
superstructure.  One  approach  to 

shipboard  HFDF  is  to  use  the 
correlation  interferometry  direction 
finding  (CIDF)  algorithm  described  by 
Saucier  and  Struckman  [6].  This 
algorithm  will  yield  accurate  bearing 
estimates  if  a  sufficiently  robust 
array  is  used.  The  algorithm  is  used 
to  compute  the  correlation  between  the 
complex  antenna  voltages  for  an 
incoming  plane  wave,  and  the  complex 
antenna  voltages  stored  in  a  data  base 
for  discrete  azimuth  angles  at  the 
same  frequency.  The  angle  of  maximum 


correlation  is  used  as  the  bearing 
estimate. 

Implementation  of  a  CIDF  system  aboard 
ship  requires  that  the  ship  be 
calibrated  to  create  the  required 
database  of  complex  receive  antenna 
voltages  for  frequencies  over  the 
range  of  interest.  This  is  not  a 
problem  as  a  ship  can  be  calibrated 
after  a  DF  system  is  installed. 
Ships,  however,  are  often  reconfigured 
as  new  systems  are  added  and  old 
systems  are  removed  or  upgraded.  The 
effect  of  associated  topside 
modifications  on  the  complex  receive 
antenna  voltages  is  generally  unknown. 
If  the  voltages  are  perturbed, 
however,  DF  system  accuracy  will  be 
degraded.  The  problem,  is  that  one 
does  not  know  if  topside  modifications 
will  cause  unacceptable  bearing 
errors . 


1.2  Problem 

At  present,  there  is  no  quantitative 
method  for  predicting  the  effect  of 
topside  reconfiguration  on  the  bearing 
accuracy  of  a  shipboard  high  frequency 
direction  finding  system  utilizing  the 
CIDF  algorithm.  One  can  recalibrate 
with  each  reconfiguration,  but  this  is 
time  consuming,  expensive,  and 
wasteful  of  ship  operating  time.  The 
objective  of  the  work  described  here 
was  to  determine  if  state— of —the— art 
hardware  and  software  would  support 
development  of  a  computer  based 
recalibration  decision  support  system. 

In  general,  one  is  interested  in  both 
HF  groundwave  and  skywave  signals,  and 
therefore  signals  which  may  be 
elliptically  polarized.  However,  when 
HF  signals  are  transmitted  from  ships 
and  submarines  at  sea,  vertically 
polarized  groundwaves  can  produce 
useful  signal  levels  at  large 
distances.  To  simplify  the  problem, 
the  scope  of  this  initial 
investigation  was  therefore  limited  to 
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vertically  polarized  ground  waves. 

1 . 3  Approach 

A  simulation  of  a  high  frequency  CIDF 
system  was  constructed  using  real  ship 
CIDF  software  to  maintain  high 
fidelity.  This  software  was  modified 
to  produce  output  files  containing 
correlation  and  bearing  estimate 
values  for  signals  tested  against  the 
database.  The  required  calibration 
database  was  created  by  computing  the 
responses  of  the  antennas  in  the  DF 
array  using  the  Numerical 
Electromagnetics  Code  (NEC) ,  version 
4.1  [7].  Wire  grid  models  of  DD963 
Spruance  Class  destroyers  with  two 
different  topside  configurations  were 
constructed  for  use  with  NEC.  One 
ship  model  had  an  Anti-Submarine 
Rocket  (ASROC)  Launcher  in  front  of 
the  deckhouse,  and  the  other  was  one 
where  the  ASROC  launcher  was  removed 
and  replaced  with  a  Vertical  Launch 
System  (VLS).  The  ASROC  launcher  was 
a  large  rectangular  box  which  was 
elevated  above  the  deck,  whereas  the 
VLS  was  recessed  into  the  deck  and  had 
a  much  lower  profile.  These  two 
configurations  represented  a  situation 
for  which  the  change  would  raise  the 
recalibration  question;  a  question 
which  at  present,  cannot  be 
satisfactorily  answered.  In  order  to 
interpret  the  large  correlation  and 
bearing  estimate  data  sets,  several 
programs  were  written  to  create  2  and 
3-dimensional  interactive  displays  of 
correlation  and  bearing  error  on  a 
Silicon  Graphics  workstation. 

NEC  4.1  was  run  on  both  a  Silicon 
Graphics  Power  Onyx  workstation  and  a 
CRAY  J90  mini  supercomputer  with  4 
processors  and  128  MWords  (1  GB)  of 
memory.  Input  and  output  files  were 
moved  over  the  School's  network 
between  these  machines  and  the  Silicon 
Graphics  workstation  on  which  the  DF 
system  simulation  was  implemented. 

Validation  was  an  important  aspect  of 
the  work  described  here.  Data  for 
validation  were  collected  in  a 
measurement  program  in  which  l/48th 
scale  brass  models  of  the  DD963  were 
used  to  measure  the  amplitude  and 
phase  responses  of  each  antenna  in  the 
direction  finding  array.  Data  were 
collected  for  both  the  ASROC  and  VLS 
configurations  at  20  frequencies  over 
the  range  88.8-1172.64  MHz  (1.85-24.43 
MHz  scaled  by  the  factor  48).  Data 


were  also  collected  for  a  single  DF 
antenna  mounted  on  a  metal  box. 

The  outdoor  measurement  facility  at 
which  the  experimental  data  were 
collected  was  located  in  San  Diego, 
CA.  The  11  foot  (3.35  meter)  long 
DD963  models  were  placed  on  a  lead 
(a=5xl06  S/m)  turntable  and  illuminated 
using  a  log  periodic  transmit  antenna. 
This  antenna  was  located  at  a  distance 
of  79  feet  (24.08  meters)  from  the 
center  of  the  turntable  and  was 
mounted  on  an  arch  which  permitted  the 
elevation  angle  to  be  varied.  A 
signal  source  and  a  vector  receiver 
were  located  in  a  control  room 
adjacent  to  the  turntable.  A  sample 
of  the  signal  from  the  source  was  fed 
to  the  receiver  for  use  as  a  phase 
reference.  Amplitude  measurements 
were  referenced  to  a  X/4  monopole  which 
was  used  for  calibration  prior  to 
recording  data.  Data  were  collected 
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The  antennas  used  in  the  DF  array  on 
the  l/48th  scale  brass  model  of  the 
DD963  were  semi-loops  with  a  230  mil 
(0.584  cm.)  radius.  These  were 
constructed  from  0.085  inch  (0.216 
cm),  outside  diameter,  semi-rigid 
coaxial  cable.  The  feed  was  a  15  mil 
(.38  mm.)  slit  cut  in  the  outer 
conductor  of  the  coax  at  the  center  of 
the  semi-loop  (highest  point  above  the 
mounting  plate)  .  At  one  end  of  the 
semi-loop,  the  coax  was  shorted  to  the 
mounting  plate,  and  at  the  other  end, 
the  coax  passed  through  the  mounting 
plate  to  provide  an  output  line. 

Accurate  computation  of  the  responses 
of  the  antennas  in  the  DF  array  was 
the  most  critical  element  of  the  work 
described  here.  Before  computing  the 
response  of  many  antennas  arrayed  on  a 
ship,  the  response  of  a  single  antenna 
mounted  on  a  metal  box,  situated  on  a 
ground  plane,  was  investigated. 

2.1  Response  of  a  Semi-loop  Antenna 
on  a  Metal  Box 

A  single  semi-loop  antenna  was  mounted 
on  one  face  of  a  metal  cube,  midway 
between  its  sides  and  8  inches  above  a 
ground  plane.  The  metal  cube  was  one 
foot  on  a  side  and  the  semi-loop  was 
contained  in  a  vertical  plane  normal 
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to  the  face  of  the  box  on  which  it  was 
mounted.  The  semi— loop  was  therefore 
responsive  to  theta  polarized  signals, 
which  includes  the  case  of  vertically 
polarized  ground  waves  which  were  of 
interest  in  this  study.  A  wire  grid 
computer  model  of  the  structure  is 
shown  in  Fig.  1,  where  it  can  be  seen 
that  the  semi-loop  was  modeled 
numerically  using  5  segments  (the 
center  segment,  #3,  is  the  feed). 
Measurements  and  computations  were 
carried  out  over  the  scaled  frequency 
range  96-1440  MHz  (2-30  MHz  unsealed) 
and  the  numerical  and  experimental 
data  were  compared.  In  each  case,  the 
excitation  was  a  theta  polarized  plane 
wave,  incident  at  an  elevation  angle 
of  10  degrees.  Fig.  2  shows  the 
amplitude  and  phase  response  at  96  MHz 
where  results  obtained  using  NEC  are 
seen  to  be  in  good  agreement  with 
experimental  results.  Note  that  the 
pattern  amplitude  is  plotted  in  dB  on 
a  relative  scale  (maximum  response  =  0 
dB).  Also  shown,  are  numerical 

results  obtained  using  the  code  PATCH 
[8].  It  should  be  noted  that 

measurements  and  computations  were 
made  for  a  range  of  elevation  angles 
with  good  agreement  observed  in  all 

caSes.  A  more  detailed  discussion  of 
this  structure,  and  the  results 
obtained  appears  in  Ref.  [9]. 

Obtaining  successful  numerical  results 
for  a  single  semi-loop  mounted  on  a 
box  gave  confidence  that  the  response 
of  the  same  antenna  could  be 
successfully  computed  when  located  on 
a  ship. 

2.2  Responses  of  the  Semi-loop 
Antennas  in  the  Shipboard  DF  Array 

The  amplitude  and  phase  responses  of 
semi-loop  antennas  in  a  shipboard  DF 
array  were  computed  at  four 

frequencies,  1.85,  6.34,  9.25,  and 

20.075  MHz,  and  compared  with 

measurement  data  at  the  scaled  values 
of  these  frequencies  (88.80,  304.32, 

444.00,  963.60  MHz).  Again,  the 

excitation  was  a  theta  polarized  plane 
wave,  but  for  the  ship,  the  elevation 
angle  was  set  at  5  degrees.  An 
elevation  angle  of  0  degrees  could  not 
be  used  for  measurements  because  the 
measurement  range  employed  a  log 
periodic  dipole  array  which  was 
mounted  on  an  arch  and  it  could  not  be 
lowered  into  the  ground.  ASROC  and 
VLS  configurations  of  the  DD963  with 
identical  24  element  DF  arrays  were 
studied.  The  DF  antennas  were 


deployed  around  the  periphery  of  the 
ship,  some  on  the  edge  of  the  deck, 
and  some  on  bulkheads,  in 
approximately  mirror  image  port  and 
starboard  locations.  Fig.  3a  shows  a 
visualization  of  the  numerical  model 
(ASROC  configuration),  and  Fig.  3b 
shows  the  l/48th  scale  brass  model 
(VLS  configuration) . 

The  numerical  models  of  the  DD963 
consisted  of  approximately  4000  wires 
and  7000  segments.  Maximum  segment 
length  was  1  meter,  or  0.1X  at  30  MHz. 
Some  areas  of  the  ship  were  meshed 
using  aim2  grid  while  other  areas 
such  as  the  sides  of  the  ship,  were 
represented  using  vertical  wires 
spaced  by  1  meter.  This  was 
considered  acceptable  since  interest 
here  was  focused  on  vertically 
polarized  signals  incident  over  the 
surface  of  the  sea.  For  convenience, 
and  to  minimize  the  number  of  segments 
in  the  numerical  model,  the  DF 
antennas  were  modeled  as  1  m2  loops  in 
exactly  the  same  locations  as  the 
semi-loop  antennas  in  the  DF  array  on 
the  l/48th  scale  brass  model.  The 
outboard  segment  of  each  loop  was  used 
as  the  feed.  Although  the  scaled  area 
of  the  square  loops  was  about  4  times 
greater  than  the  area  of  the  semi¬ 
loops,  both  were  electrically  small, 
so  the  difference  was  of  no 
consequence.  There  were  several 
areas  along  the  sides  of  the  ship 
(Fig.  3a)  between  the  deck  and  the 
waterline,  however,  where  there  were 
openings  as  large  as  15  meters.  As 
will  be  shown  later,  these  areas 
appear  to  have  caused  error  in  the 
numerical  results  at  the  higher 
frequencies . 

Numerical  and  experimental  values  of 
amplitude  and  phase  were  generated  for 
each  of  the  24  antennas  in  the  DF 
array  at  the  four  frequencies 
mentioned  earlier.  This  resulted  in 
far  more  data  than  can  be  presented 
here,  so  only  representative  results 
will  be  shown.  The  discussion  which 
follows  refers  to  specific  antennas  in 
the  DF  array.  The  antennas  are 
numbered  1  through  12  going  from  bow 
to  stern  and  are  designated  as  port 
(P)  or  starboard  (S) . 

To  facilitate  comparison  of  data,  the 
results  have  been  displayed  for  port 
and  starboard  antenna  pairs,  with 
results  for  both  ASROC  and  VLS 
configurations  shown  on  the  same  plot. 
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Pattern  amplitudes  have  been 
normalized  (maximum  response  =  0  dB) 
to  make  it  easier  to  compare 
differences  between  the 
configurations.  The  phase  reference 
is  arbitrary  both  numerically  and 
experimentally,  and  no  attempt  was 
made  to  match  the  references.  The 
numerical  model  of  the  DD963  is 
oriented  along  the  x  axis  in  a  polar 
spherical  coordinate  system.  Azimuth 
is  determined  by  the  angle  <pr  which  is 
measured  in  the  x-y  plane,  in  the 
counter-clockwise  direction  from  the  x 
axis,  in  accordance  with  mathematical 
convention.  For  measured  data,  angle 
is  measured  clockwise  from  the  bow  in 
accordance  with  nautical  convention. 
In  each  case,  however,  0  degrees 
corresponds  to  the  bow  of  the  ship, 
and  the  antenna  patterns  relative  to 
the  bow  of  the  ship  are  correctly 
depicted,  with  the  port  antenna  on  the 
left  and  the  starboard  antenna  on  the 
right.  Thus  numerical  and 
experimental  patterns  can  be  easily 
compared  and  their  spatial  responses 
are  correctly  depicted.  In  the 
following  discussion,  the  term 
"azimuth"  means  relative  azimuth  in 
the  nautical  sense. 

Antennas  P-3  and  S-3  are  located  on 
the  front  of  the  deckhouse,  just 
behind  the  ASROC  launcher  (Fig.  3a)  . 
Since  these  antennas  are  closest  to 
the  ASROC  launcher,  it  was  anticipated 
that  they  would  be  most  affected  by 
its  removal.  The  surfaces  on  which 
these  two  antennas  are  mounted  are 
angled  outboard  by  about  20  degrees. 
Thus,  these  antennas  "look"  20  degrees 
to  port  and  starboard  of  the  bow, 
respectively. 

Figs.  4a  and  4b  show  the  numerical  and 
experimental  patterns  for  antennas  P-3 
and  S-3  at  1.85  MHz.  Results  for  the 
ASROC  configuration  appear  as  a  solid 
curve  and  results  for  the  VLS 
configuration  appear  as  a  dashed 
curve.  Because  the  curves  are 
reproduced  here  in  black  and  white 
rather  than  color,  the  two  patterns 
are  indistinguishable  in  some  cases. 
In  this  case,  the  patterns  for  the  two 
configurations  are  the  same, 
indicating  that  the  removal  of  the 
ASROC  launcher  had  little  affect  on 
the  patterns  of  these  antennas  at  this 
frequency.  Both  numerical  and 
experimental  results  show  the  patterns 
are  approximately  omnidirectional  with 
a  maximum  response  about  20  degrees 


off  the  bow.  The  agreement  is  quite 
good. 

Figs.  5a  and  5b  show  the  numerical  and 
experimental  phases  for  antennas  P-3 
and  S-3  at  1.85  MHz.  Again,  there  is 
essentially  no  difference  between  the 
results  for  the  ASROC  and  VLS 
configurations.  Phase  is  plotted 
within  a  360  degree  range,  so  where 
phase  switching  occurs,  one  must 
visualize  an  n27r  (n=integer) 
translation  which  would  create  a 
continuous  curve.  Also,  when 
comparing  numerical  and  experimental 
results,  the  fact  that  the  phase 
references  are  arbitrary,  and 
different,  must  be  taken  into  account. 
Thus,  numerical  and  experimental 
curves  must  be  overlaid  and  shifted 
vertically  to  achieve  alignment  so 
they  can  be  compared.  This  can  be 
done  easily  at  low  frequencies,  as  in 
the  case  of  Figs.  5,  but  visual 
comparison  is  difficult  at  higher 
frequencies  where  phase  varies  rapidly 
and  there  are  many  360  degree 
discontinuities . 

The  numerical  and  experimental  phase 
responses  shown  in  Figs.  5a  and  5b  are 
in  reasonable  agreement.  Numerically, 
the  total  phase  variation  is  about  175 
degrees.  Experimentally,  it  is  about 
135  degrees.  Most  of  the  phase  change 
which  occurs  for  a  DF  antenna  is  due 
to  the  offset  of  the  antenna  from  the 
center  of  rotation,  which  is  the 
center  of  the  ship.  Experimentally, 
for  example,  when  the  ship  is  rotated 
on  the  turntable,  a  given  DF  antenna 
will  move  toward  or  away  from  the 
source  as  the  ship  rotates.  This  will 
cause  phase  to  advance  or  retard 
relative  to  a  reference  of  fixed 
phase.  Antennas  P-3  and  S-3,  for 
example,  are  located  a  full  scale 
distance  of  39.56  meters  from  the 
center  of  rotation.  Thus,  the 
distance  of  these  antennas  from  the 
source  changes  by  79.1  meters  or  0.49X 
at  1.85  MHz.  This  change  in  position 
corresponds  to  a  phase  change  of  175 
degrees,  which  is  in  good  agreement 
with  the  numerical  result. 

The  experimental  results  for  P-3  and 
S-3  are  offset  from  one  another  by  180 
degrees.  The  phase  curve  for  P-3 
starts  at  about  155  degrees,  while 
that  for  S-3  starts  at  about  335 
degrees.  These  results  indicate  that 
on  the  brass  model,  antenna  P-3  was 
mounted  upside  down  (rotated  180 
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degrees)  with  respect  to  S-3.  This  is 
interesting,  but  has  no  effect  on  the 
correlation  and  bearing  error  results 
presented  later.  It  would  have  an 
impact,  however,  if  numerical  and 
experimental  results  were  cross- 
correlated. 

Figs,  6a  and  6b  show  the  numerical  and 
experimental  patterns  for  DF  antennas 
P-3  and  S-3  at  9.25  MHz.  At  this 
frequency,  the  numerical  results  show 
the  port  antenna  looking  about  20 
degrees  to  starboard  and  the  starboard 
antenna  looking  about  15  degrees  to 
port  (the  ship  is  not  symmetric  about 
the  keel).  The  experimental  patterns 
show  some  scalloping  in  these 
directions  and  maxima  to  either  side. 
The  numerical  patterns  show  a  lower 
response  in  the  stern  sector  than  the 
experimental  patterns.  The 
experimental  results  show  more 
difference  between  the  ASROC  and  VLS 
configurations.  The  agreement  between 
numerical  and  experimental  results  at 
this  frequency  is  still  reasonable, 
although  less  satisfying  than  at  1.85 
MHz,  probably  due  to  deterioriating 
fidelity  of  the  numerical  model. 

As  a  last  example,  Figs.  7  a  and  7b 
show  the  patterns  of  DF  antennas  P-12 
and  S-12  at  20.075  MHz.  These 
antennas  are  mounted  on  the  stern, 
looking  aft.  At  this  frequency,  the 
numerical  and  experimental  patterns 
for  some  antennas  are  still  in 
reasonable  agreement,  but  for  others 
there  is  significant  error.  For  P-12 
and  S-12,  the  experimental  results 
show  a  maximum  response  at  180  degrees 
relative  azimuth,  as  one  would  expect, 
and  response  in  the  bow  sector  is  down 
15-30  dB.  The  numerical  results  are 
quite  inaccurate  for  these  antennas. 
The  reason  for  this,  it  is  suspected, 
is  that  the  several  large  openings  in 
the  side  of  the  hull  (d=\;  see  Fig.  3a) 
resulted  in  excitation  of  a  field 
inside  the  hull  at  20.075  MHz. 

It  is  not  possible  to  present  all  of 
the  numerical  and  experimental  data 
which  were  generated  in  the  course  of 
this  study,  but  those  data  which  have 
been  presented  above  are 
representative.  Related  results  can 
be  found  in  Ref.  [10].  The  conclusion 
reached,  based  on  inspection  of  all 
numerical  and  experimental  results,  is 
that  the  numerical  models  of  the  DD963 
produce  good  results  up  to  about  9 
MHz,  but  need  improvement  for  use  at 


higher  frequencies.  The  addition  of 
vertical  wires  to  close  the  openings 
in  the  sides  of  the  hull  would 
probably  result  in  significant 
improvement  and  will  be  investigated. 

3.  Correlation  and  Bearing  Error 

Numerical  and  experimental  DF  antenna 
amplitude  and  phase  data  were  used  in 
the  computer  simulation  of  the  DD963 
DF  system  to  study  the  impact  of 
topside  configuration  changes  on 
system  performance.  As  mentioned 
earlier,  the  specific  configuration 
change  examined  in  the  work  described 
here  was  the  removal  of  an  ASROC 
launcher  (located  just  in  front  of  the 
deckhouse)  and  its  replacement  with  a 
VLS. 

The  DF  antenna  responses  (amplitude 
and  phase  vs.  azimuth)  for  the  ASROC 
configuration  constitute  a  database 
for  that  configuration.  The  DF 
antenna  responses  for  the  VLS 
configuration  constitute  a  database 
for  that  configuration,  and  also 
represent  the  responses  which  would  be 
measured  if  a  signal  were  incident  on 
the  ship  in  that  configuration.  Thus, 
if  a  signal  vector  (the  set  of  complex 
DF  antenna  responses  for  a  given 
azimuth)  for  the  VLS  configuration  is 
cross-correlated  with  signal  vectors 
for  the  ASROC  configuration  over  all 
azimuths,  the  result  is  a  cross 
correlation  curve,  the  peak  of  which 
may  be  used  to  derive  an  estimate  of 
the  signal  angle  of  arrival.  This 
procedure  simulates  the  effect  of 
calibrating  the  ship  in  the  ASROC 
configuration,  then  changing 
configuration  to  VLS,  and  using  the 
ASROC  database  to  DF  signals  received 
in  the  VLS  configuration. 

Numerical  and  experimental  DF  antenna 
response  data  were  used  to  create 
databases  for  the  ASROC  and  VLS 
configurations  of  the  DD963  at  four 
frequencies;  1.85,  6.34,  9.25,  and 
20.075  MHz.  The  VLS  databases  were 
cross-correlated  with  the  ASROC 
databases  to  determine  the  magnitude 
of  the  complex  cross-correlation 
coefficient,  derive  a  bearing 
estimate,  and  compute  bearing  error. 
These  data  were  then  used  to  generate 
several  displays  which  permitted  the 
results  to  be  interpreted. 
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3 . 1  Cross-correlation  Surfaces 

Figs*  8a  and  8b  show  the  numerical  and 
experimental  surfaces  of  cross¬ 
correlation  for  1.85  MHz.  The 
vertical  (z)  axis  is  the  magnitude  of 
the  complex  correlation  coefficient 
( 0< | R | <1 ) ,  the  horizontal  (x)  axis  is 
the  difference  between  the  azimuth 
angle,  <p,  of  the  VLS  signal  and  the 
azimuth  angle,  <plf  of  the  ASROC 
database  signal  (-18O°<(0-0t) <+180°  ) , 
and  the  axis  receeding  into  the  page 
(y  axis)  is  the  azimuth  angle  of 
arrival,  <£,  of  the  VLS  signal 
(O°<0<36O°).  If  the  DF  system  is 
performing  properly,  the  surface  will 
have  low  sidelobes,  and  a  central 
ridge  near  (0~0t)=O  with  |r|~1.  A 
configuration  change  which  causes 
serious  bearing  errors  will  result  in 
a  surface  with  maximum  correlation  at 
other  locations  (0-0t)^O,  usually  with 
| R | <1 . 

Although  only  2D  pictures  can  be 
presented  here,  in  the  actual 
simulation,  the  computer  display  of 
the  cross-correlation  surface  is 
interactive  and  can  be  rotated  about 
any  axis  to  examine  its  properties. 
Figs.  8a  and  8b  show  that  simulation 
results  using  numerical  and 
experimental  DF  antenna  response  data 
for  the  ASROC  and  VLS  configurations 
agree  quite  well  at  1.85  MHz. 

3.2  Cross-correlation  vs.  Azimuth 

A  cut  through  the  cross-correlation 
surface,  perpendicular  to  the  <p  axis, 
yields  a  cross-correlation  curve  for  a 
fixed  signal  angle  of  arrival,  <p.  A 
polar  plot  of  this  curve  gives  the 
relative  response  of  the  DF  array  vs. 
azimuth,  in  a  format  similar  to  that 
of  a  normal  antenna  pattern.  Figs.  9a 
and  9b,  show  cuts  through  the 
numerical  and  experimental  surfaces 
(Figs.  8a  and  8b)  for  a  signal  arrival 
angle  of  232  degrees,  as  indicated  by 
the  radial  cursor.  The  numerical  and 
experimental  results  both  show  peak 
correlation  near  232  degrees, 
indicating  little  bearing  error  for 
this  angle  of  arrival. 

Again,  only  a  static  2D  result  can  be 
presented  here,  but  in  the  actual 
simulation,  the  radial  cursor,  which 
indicates  the  angle  of  arrival  on  the 
computer  display,  can  be  moved  to  any 
azimuth  using  a  mouse.  As  the  cursor 
moves  on  the  computer  screen,  the 


proper  cross-correlation  curve  is 
displayed.  Thus,  a  dynamic  picture  of 
the  array  response  vs.  azimuth  can  be 
quickly  and  easily  obtained,  and 
sectors  where  bearing  errors  occur  can 
be  easily  identified. 

3.3  Bearing  Error  vs.  Azimuth 

Figs.  10a  and  10b  show  bearing  error 
vs.  azimuth  at  1.85  MHz  as  determined 
from  the  simulation  using  numerical 
and  experimental  DF  antenna  responses. 
The  curves  both  show  a  worst  case 
bearing  error  of  approximately  ±2.5°. 
In  particular,  at  232°,  both  curves 
show  a  small  negative  bearing  error. 
This  negative  error  (bearing  estimate 
less  than  true  bearing)  is  also 
evident  in  Figs.  9a  and  9b  where  the 
peak  correlation  can  be  seen  to  occur 
at  a  bearing  slightly  less  than  that 
of  the  cursor  which  indicates  the  true 
bearing. 

The  simulation  can  also  be  used  to 
create  a  display  of  the  3  dimensional 
surface  of  bearing  error  vs.  azimuth 
and  frequency  from  curves  of  the  type 
shown  in  Figs.  10.  This  surface  makes 
it  possible  to  easily  identify  areas 
of  high  bearing  error  over  the  entire 
domain  of  azimuth  and  frequency. 

3 • 4  RMS  Bearing  Error 

The  last  display  which  can  be  created 
using  the  simulation  is  one  of  RMS 
bearing  error  vs.  frequency.  In  this 
display,  the  root  mean  square  bearing 
error  is  computed  over  all  azimuth 
angles  at  each  frequency  of  interest. 
This  provides  an  integrated  measure  of 
DF  system  performance.  Table  1 


Freq. 

(MHz) 

1.85 

6.34 

9.25 

20.077 

Num. 

1.11° 

0.86° 

2.80° 

7.38° 

Expt . 

1.20° 

0.39° 

1.51° 

6.20° 

Table  1.  RMS  Bearing  Error  Predicted 
by  the  Simulation. 

shows  the  results  obtained  using 
numerical  and  experimental  DF  antenna 
data  at  the  four  frequencies 
investigated  in  this  study. 

As  indicated  earlier,  the  accuracy  of 
the  numerical  values  of  DF  antenna 
response  degraded  with  increasing 
frequency  due  to  (most  probably 
correctable)  deficiencies  in  the  wire 
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grid  model  of  the  DD963.  This 
conclusion  was  reached  by  comparing 
numerical  and  experimental  results  and 
noting  that  the  experimental  antenna 
patterns  appeared  as  one  would  expect 
based  on  frequency  and  location  while 
the  numerical  patterns  at  the  higher 
frequencies  did  not.  That  is,  there 
is  no  reason  to  suspect  that  any  of 
the  experimental  data  is  incorrect. 
Thus,  simulation  results  (correlation, 
bearing  error)  based  on  numerical  and 
experimental  data  were  also  less 
consistent  at  the  higher  frequencies. 
This  was  expected  as  the  simulation 
output  depended  on  the  DF  antenna 
data.  However,  the  data  in  Table  1 
shows  that  even  with  this  degradation, 
the  RMS  bearing  error  results  track 
quite  well  since  this  is  a  performance 
measure  which  is  integrated  over 
azimuth.  One  can  therefore  tenatively 
conclude  that  the  computer  simulation 
will  correctly  predict  when  bearing 
error  will  be  caused  by  a  topside 
change  of  configuration  even  at 
frequencies  where  there  is 
considerable  error  in  the  computation 
of  some  DF  antenna  responses. 

4 •  Conclusions 

This  paper  has  presented  the  results 
of  a  study  to  determine  if  computer 
simulation  could  be  used  to  build  a 
decision  support  system  to  determine 
when  topside  changes  in  configuration 
might  require  recalibration  of  a 
ship's  CIDF  based  high  frequency 
direction  finding  system.  A  computer 
simulation  of  the  DF  system  was 
constructed,  DF  antenna  responses 
required  to  create  the  system  database 
were  computed,  and  measurements  were 
made  for  validation  purposes. 

The  most  critical  aspect  of  this  work 
was  the  use  of  NEC  4.1  to  compute  the 
DF  antenna  responses.  Satisfying 
results  were  obtained  at  the  lower 
frequencies,  but  computed  antenna 
patterns  became  less  accurate  with 
increasing  frequency.  This  appears  to 
be  due  to  ship  numerical  modeling 
deficiencies  which  are  believed  to  be 
correctable.  Thus,  the  approach  seems 
quite  promising. 

The  need  for  additional  work  on 
several  enabling  technologies  quickly 
became  clear  during  the  course  of  this 
investigation.  First,  methods  which 
have  previously  been  used  to  create 
numerical  models  of  complex  structures 


are  laborious  and  time  consuming.  For 
a  computer  simulation  such  as  that 
described  in  this  paper  to  be 
realistically  implemented,  it  must  be 
possible  to  create  and  modify 
numerical  models  quickly  and  easily. 
Ideally,  one  would  like  to  go 
seamlessly  from  a  ship  CAD  file  to  the 
input  file  for  a  CEM  code.  Presently, 
one  must  use  software  to  mesh  a 
structure  of  interest,  fix  the  result 
manually,  and  then  translate  the 
output  to  obtain  an  input  file  for  the 
CEM  code.  When  the  code  is  run, 
errors  will  likely  occur  and  again, 
manual  work  will  be  required  to  fix 
these.  More  work  is  needed  here. 

Good  pre-processors  and  post¬ 
processors  for  NEC  (or  any  other  code) 
are  also  needed.  There  is  a  need  to 
be  able  to  quickly  and  easily  find 
wires  and  segments  identified  in  CEM 
code  error  and  warning  messages  and  to 
eliminate  the  conditions  causing  the 
problems.  Model  visualization  and 
file  editing  capabilities  are  being 
developed  but  have  not  reached  the 
point  where  this  is  a  straight  forward 
task  for  a  large  file.  Lastly, 
additional  work  is  needed  on  post 
processor  display  capabilities.  The 
processor  used  for  this  work,  for 
example,  did  not  permit  display  of 
receive  antenna  amplitude  and  phase 
data  from  the  NEC  output  file,  and  it 
was  therefore  necessary  to  develop 
separate  software  for  this  purpose. 

Lastly,  difficulties  were  encountered 
using  NEC  to  determine  the  receiving 
patterns  of  multiple  antennas.  For 
this,  one  needs  the  antenna  feedpoint 
currents.  To  print  only  the  24 
antenna  feedpoint  currents,  it  was 
necessary  to  place  the  feedpoint 
segments  at  the  beginning  of  the  file, 
so  they  appeared  in  sequence  and  the 
NEC  print  command  (PT)  could  be  used 
to  obtain  all  24  currents  in  one  run. 
This  required  breaking  up  the  input 
file  which  is  undesirable.  It  would 
be  much  more  desirable  to  be  able  to 
reuse  the  PT  command  to  force  printing 
of  the  current  for  each  segment 
desired  in  a  given  run. 
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Figure  2.  Numerical  and  experimental  amplitude  and  phase  responses  for  semi¬ 
loop  antenna  on  a  metal  box.  Frequency  is  96  MHz  (2  MHz  scaled  X48).  (  = 

measured,  —  =  PATCH,  .  .  =  NEC) 


DD963:  NEC4  Antenna  RX  Patterns 
Frequency  =  1.85  MHz 

Eta  =  0  Deg  Theta  =  85  Deg  Elevation  =  5  Deg 
LEGENDS 


ASROC  VLS  ASROC  VLS 

Norm.  Factor  0.02658  0.02732  Norm.  Factor  0.02725  0.02799 

Data  Files...  asrocf185rx.dat  vlsf185ixdat  Data  Files...  asrocf185rx.dat  vlsf185rxdat 

Segment  #6  6  Segment  #5  5 

Figure  4a,  Numerical  patterns  of  DF  antennas  P-3  and  S-3  for  ASROC  and  VLS 
configurations  of  DD963*  Elevation  angle  is  5  degrees,  frequency  is  1*85  MHz* 


DD963  Brass  Model:  Measured  Antenna  RX  Patterns 
Frequency:  1.85  MHz 
Elevation:  5  Deg 
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Figure  4b.  Experimental  patterns  of  DF  antennas  P-3  and  S-3  for  ASROC  and  VLS 
configurations  of  DD963.  Elevation  angle  is  5  degrees,  scaled  freguencv  is 
1.85  MHz. 
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Figure  5a.  Numerical  phase  of  DF  antennas  P-3  and  S-3  for  ASROC  and  VLS 
configurations  of  DD963.  Elevation  angle  is  5  degrees,  frequency  is  1.85  MHz. 
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Figure  5b.  Experimental  phase  of  antennas  P-3  and  S-3  for  ASROC  and  VLS 
configurations  of  DD963.  Elevation  angle  is  5  degrees,  scaled  frequency  is 
1.85  MHz. 
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Figure  6a,  Numerical  patterns  of  antennas  P-3  and  S-3  for  ASROC  and  VLS 
configurations  of  DD963*  Elevation  angle  is  5  degrees,  frequency  is  9,25  MHz 
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Figure  6b.  Experimental  patterns  of  antennas  P-3  and  S-3  for  ASROC  and  VLS 
conf  icju^stions  of  DD963 .  Elevation  angle  is  5  degrees,  scaled  frecruencv  is 
9.25  MHz. 
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Figure  7a.  Numerical  patterns  of  DF  antennas  P-12  and  S-12  for  ASROC  and  VLS 
configurations  of  DD963.  Elevation  angle  is  5  degrees,  frequency  is  20.075 
MHz. 
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Figure  7b.  Experimental  patterns  for  DF  antennas  P-12  and  S-12  for  ASROC  and 
VLS  configurations  of  DD963.  Elevation  angle  is  5  degrees,  scaled  frequency 
is  20.075  MHz. 


Figure  8.  Cross-correlation  surface,  VLS  vs.  ASROC  configuration  of  DD963. 
(a)  Numerical  result  (top),  1.85  MHz,  (b)  Experimental  result  (bottom),  1.85 
MHz  scaled. 
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Abstract 

A  method  of  moments  (MoM)  code  has  been 
developed  to  compute  the  scattering  from  a  pla¬ 
nar  or  cylindrically  conformed  slot  array  antenna. 

By  hybridizing  the  MoM  with  the  shooting  and 
bouncing  ray  (SBR)  method,  the  scattering  from 
a  large,  complex  target  with  a  slot  array  antenna 
can  be  computed.  The  scattering  problem  can 
be  decomposed  using  the  field  equivalence  prin¬ 
ciple  such  that  the  MoM  is  employed  to  model 
the  slot  array  while  the  SBR  method  is  used  to 
compute  the  scattering  from  the  large,  complex 
target.  Sample  results  show  the  utility  of  the 
method  and  the  need  to  include  slot  array  scat¬ 
tering  when  computing  the  RCS  of  a  complex 
target. 

1  Introduction 

The  presence  of  a  slotted  waveguide  array  antenna  on  a 
radar  target  may  have  a  significant  contribution  to  the 
overall  radar  cross-section  (RCS)  of  the  target.  There¬ 
fore,  the  computation  of  the  RCS  should  include  the 
scattering  from  the  slot  array.  Recently,  a  method  of  mo¬ 
ments  (MoM)  procedure  has  been  introduced  to  compute 
the  scattering  from  a  cylindrically  conformal  slotted- 
waveguide  array  antenna  [1,2].  However,  this  procedure 
does  not  take  into  account  the  geometry  in  which  the 
slot  array  is  located.  If  the  slot  array  is  located  in  a 
complex,  three-dimensional  (3-D)  geometry,  the  MoM 
cannot  efficiently  account  for  the  effect  of  the  geome¬ 
try.  A  more  efficient  method  to  compute  the  scattering 
from  a  large,  3-D  body  is  the  high  frequency  shooting- 
and-bouncing-ray  (SBR)  method.  However,  this  method 
cannot  accurately  account  for  the  slots,  each  of  which  is 
typically  smaller  than  an  electromagnetic  wavelength  in 
size.  In  this  paper,  the  MoM  computation  of  the  scatter¬ 
ing  from  a  slot  array  is  hybridized  with  the  SBR  method 
to  compute  the  electromagnetic  scattering  from  a  large, 
3-D  target  which  includes  a  slot  array  antenna. 

The  basis  of  the  hybrid  method  is  the  field  equiva¬ 
lence  principle,  which  allows  the  scattering  geometry  to 


be  decomposed  into  separate  regions.  The  MoM  is  ap¬ 
plied  to  the  slotted  waveguides,  while  the  SBR  method 
is  applied  to  the  region  outside  the  waveguides,  which 
includes  the  complex,  3-D  target.  By  using  the  hybrid 
method,  the  scattering  from  a  large,  3-D  target,  which 
includes  a  slotted-waveguide  array  antenna,  can  be  effi¬ 
ciently  and  accurately  computed. 

The  remainder  of  this  paper  is  divided  into  four  sec¬ 
tions.  Section  2  describes  the  formulation  of  the  prob¬ 
lem,  including  the  use  of  the  MoM,  the  use  of  the  SBR 
method,  and  techniques  to  decouple  the  computations 
of  the  two  methods.  Section  3  describes  briefly  how  the 
method  has  been  tested,  and  Section  4  gives  some  nu¬ 
merical  results  which  show  the  capability  of  the  method. 
The  results  in  Section  4  also  demonstrate  the  need  to  in¬ 
clude  the  slot  array  in  scattering  computations.  Finally, 
Section  5  gives  a  brief  conclusion. 


2  Formulation 

Consider  the  example  target  shown  in  Figure  la.  The 
target  is  complex  and  3-D,  and  it  includes  a  slotted 
waveguide  array  antenna  on  its  surface.  The  slotted 
waveguide  array  antenna  may  be  planar,  or  it  may  con¬ 
form  to  the  surface  of  a  cylinder.  The  first  step  to  com¬ 
pute  the  scattering  from  this  target  is  to  analyze  the 
slotted  waveguides  using  the  MoM.  Then,  the  scattering 
from  the  target  with  the  slot  apertures  covered  by  per¬ 
fect  electric  conductor  (PEC)  is  computed  using  the  SBR 
method.  During  the  SBR  calculation,  the  incident  field 
on  the  slot  array  antenna  is  computed  and  stored.  This 
incident  field  is  combined  with  the  MoM  analysis  to  find 
an  equivalent  magnetic  current  on  the  outer  aperture  of 
each  slot.  Finally,  the  radiation  of  these  equivalent  mag¬ 
netic  currents  in  the  presence  of  the  complex,  3-D  target 
is  computed  using  the  reciprocity  theorem.  This  result  is 
added  to  the  previously  computed  SBR  scattering  result. 


43 


— 

I  IP#;: 


Region  I 
a  (outside) 


Region  HI 
(waveguide) 


Region  II 
(slot) 


(b)  Local  slot  geometry 


Figure  1:  Example  of  a  complex,  3-D  target  with  a  slot  array  antenna  and  the  local  slot  geometry. 


2.1  Use  of  MoM 

The  first  step  in  the  formulation  of  the  problem  is  to 
analyze  the  slotted  waveguides  using  the  MoM.  There 
are  two  main  steps  in  the  application  of  the  MoM.  First, 
the  problem  must  be  described  in  terms  of  an  integral 
equation.  Then,  the  integral  equation  is  discretized  to 
find  a  numerical  solution.  The  steps  are  outlined  here, 
and  more  detail  is  given  in  [1,2]. 

To  derive  the  integral  equation,  the  apertures  of  each 
slot  are  first  covered  with  PEC,  and  equivalent  magnetic 
currents  over  each  aperture  are  introduced.  Figure  lb 
depicts  the  situation  for  the  ith  slot.  The  region  outside 
of  the  antenna  is  denoted  Region  I,  the  region  inside  the 
slot  is  Region  II,  and  the  region  outside  of  the  slot  but 
inside  the  waveguide  is  Region  III.  An  equivalent  mag¬ 
netic  current  M?  is  introduced  on  the  inside  of  the  outer 
slot  aperture  (between  Regions  I  and  II),  and  the  equiv¬ 
alent  current  is  introduced  on  the  waveguide  side 
of  the  inner  aperture  (between  Regions  II  and  III).  Be¬ 
cause  the  electric  field  must  be  continuous  across  each 
aperture,  — M f  must  be  introduced  on  the  outside  of  the 
outer  aperture,  and  — must  be  introduced  on  the  slot 
side  of  the  inner  aperture.  Note  that  when  the  analysis 
is  completed,  — are  the  currents  that  radiate  in  the 
presence  of  the  complex,  3-D  body  as  discussed  above. 

To  derive  the  integral  equation,  the  continuity  of 
the  tangential  magnetic  field  across  each  aperture  is  en¬ 
forced.  Denoting  the  tangential  magnetic  field  in  Region 
III  on  the  ith  slot  aperture  due  to  the  magnetic  current 
on  the  jth  aperture  as  Hj^Mj),  the  following  must  hold 
on  each  inner  aperture: 

E  +  HS(M?)  -  H?f(M f)  =  0.  (1) 

3 

Further,  denoting  the  tangential  incident  field  on  the  ith 


slot  aperture  as  H^fR, 

£H-'(Mi)  +  HS(M?)  -  HS(M*6)  =  H*?R 
i  (2) 

must  hold  on  each  outer  slot  aperture.  Note  that  the  in¬ 
cident  fields  are  calculated  using  the  SBR  method,  and 
the  magnetic  field  due  to  a  magnetic  current  is  found 
from 

H“(M)  =  JJ  da{r,r')  M(r’)dS’  (3) 

5 

where  a  is  I,  II,  or  III,  depending  on  the  region  of  in¬ 
terest,  G  (r,  r')  is  the  magnetic-source-magnetic-field 
dyadic  Green’s  function  in  the  appropriate  region,  and  r 
corresponds  to  the  point  at  which  the  magnetic  field  is 
to  be  evaluated.  Combining  Equations  1,  2,  and  3  gives 
an  integral  equation  for  the  magnetic  currents. 

The  second  main  step  in  application  of  the  MoM  is 
to  discretize  the  integral  equation  to  find  a  numerical  so¬ 
lution  for  the  currents.  To  accomplish  this  step,  the  cur¬ 
rents  are  expanded  in  terms  of  sinusoidal  basis  functions. 
Defining  £  to  be  the  direction  parallel  to  the  lengths  of 
the  slots  and  using  a  local  coordinate  system  in  which 
£j  =  0  at  one  end  of  the  jth  slot,  the  current  on  the  jth 
slot  aperture  is  expanded  as 

M;=«XXsi”(f;D  H> 

where  N  is  the  number  of  terms  in  the  expansion,  and 
/?  represents  a  for  the  current  on  the  outer  aperture  or  b 
for  the  current  on  the  inner  aperture.  Equation  4  is  valid 
for  points  on  the  jth  slot  aperture;  for  points  outside  of 
the  aperture,  the  expansion  is  defined  to  be  zero.  As¬ 
suming  the  width  of  a  slot  is  much  less  than  its  length, 
the  £  component  of  the  current  is  the  only  component  of 
interest. 
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Substituting  the  expansion  given  in  Equation  4  into 
the  integral  equation  allows  the  integral  equation  to  be 
converted  to  a  matrix  equation  which  can  be  solved  nu¬ 
merically.  For  more  details  about  solving  the  integral 
equation,  the  reader  is  referred  to  [1,2].  However,  one 
important  step  that  should  be  mentioned  here  is  the 
derivation  of  the  dyadic  Green’s  functions  for  the  var¬ 
ious  regions.  The  Green’s  functions  given  in  [1,2]  for 
Regions  II  and  III  are  applicable  to  the  present  prob¬ 
lem.  For  Region  I,  the  dyadic  Green’s  function  can  be 
written  as 

_ t  —  cvl  ■ 1  diff 

G  (r,  r')  =  G  (r,r')  +  G  (r,r').  (5) 

The  Green’s  function  given  in  [1,2]  for  the  exterior  re- 

— cy  1  ~  diff 

gion  corresponds  to  G  (r,rr),  and  G  (i*,!*7)  is  a  per¬ 
turbation  term  due  to  diffraction  and  reflection  by  the 
complex  target  in  which  the  slot  array  is  embedded.  Ne¬ 
glecting  (r,  r')  neglects  fields  which  are  scattered  by 
the  slots,  diffracted  or  reflected  by  the  large  body  back 
to  the  slots,  and  scattered  by  the  slots  again  [3].  These 
fields  are  usually  an  insignificant  part  of  the  scattering, 
and  this  term  is  neglected  in  the  computations.  Thus, 
the  Green’s  function  given  in  [1,2]  for  Region  I  is  used 
for  the  present  problem. 

2.2  Use  of  SBR 

As  previously  mentioned,  the  MoM  is  used  to  analyze 
the  slot  array  antenna  while  the  SBR  method  is  used 
for  the  remainder  of  the  problem.  Thus,  there  are  three 
main  tasks  to  be  accomplished  by  the  SBR  method:  to 
compute  the  scattering  from  the  complex,  3-D  target,  to 
compute  the  incident  magnetic  fields  on  the  slot  aper¬ 
tures,  and  to  compute  the  radiation  of  the  equivalent 
currents  on  the  slot  apertures  in  the  presence  of  the  com¬ 
plex,  3-D  target.  In  all  of  these  cases,  the  slot  apertures 
are  covered  with  PEC. 

The  details  of  the  SBR  method  are  discussed  in  [3- 
6].  For  the  present  problem,  the  SBR  procedure  is  used 
to  compute  the  scattering  from  the  complex,  3-D  tar¬ 
get  with  the  slot  apertures  closed  by  PEC.  The  SBR 
procedure  is  implemented  using  the  XPATCH  software 
package  [4,5]. 

The  incident  magnetic  field  on  the  slot  apertures  is 
computed  using  SBR  at  the  same  time  the  scattering 
from  the  complex,  3-D  target  is  computed.  While  trac¬ 
ing  the  rays  to  find  the  scattering,  some  rays  will  hit  on 
or  near  a  slot  aperture.  The  field  contributions  from  each 
of  these  rays  are  combined  with  appropriate  phase  shifts 
to  find  the  incident  magnetic  field  on  each  slot  aperture. 
The  incident  magnetic  fields  on  the  slot  apertures  are 
used  by  the  MoM  to  compute  the  equivalent  magnetic 
currents  on  the  apertures. 


The  remaining  step  in  the  problem  is  to  compute 
the  radiation  of  the  magnetic  currents  in  the  presence  of 
the  large  body.  The  SBR  method  together  with  the  reci¬ 
procity  theorem  is  employed  for  this  task  [3,6].  Consider 
an  infinitesimal  dipole  placed  at  the  scattering  observa¬ 
tion  point.  If  the  target  containing  the  slot  array  is  in  the 
far  field  of  the  dipole,  the  dipole  launches  a  plane  wave 
toward  this  target.  Recall  that  for  the  SBR  method,  the 
grid  of  rays  launched  toward  the  target  corresponds  to 
a  plane  wave.  Note  also  that  the  reciprocity  theorem 
states 

Eslot  •  JdV  =  j  j  H®br  •  MadS  (6) 
v  s 

where  H®BR  is  the  incident  field  on  the  slot  apertures  due 
to  the  dipole  at  the  scattering  observation  point,  Ma  is 
the  current  on  the  outer  slot  apertures,  which  is  found 
using  the  MoM,  E^0^  is  the  radiation  due  to  ,  and 
J  is  the  dipole  current.  Thus,  if  the  dipole  current  (J)  is 
appropriately  chosen  and  mono-static  scattering  is  being 
computed,  all  components  to  find  Eslot  using  reciprocity 
are  computed  already.  If  bi-static  scattering  results  are 
desired,  H®BR  resulting  from  a  dipole  at  the  scattering 
observation  point  must  be  computed  first,  then  Eslot  can 
be  computed. 

2.3  Decoupling  the  MoM  from  the  SBR 
Method 

As  they  are  presented  in  Section  2.1,  the  MoM  compu¬ 
tations  are  coupled  to  the  SBR  method  computations. 
This  is  due  to  the  fact  that  the  incident  magnetic  field 
on  the  slot  apertures,  which  is  computed  using  the  SBR 
method,  is  required  for  the  MoM  computations.  To  avoid 
having  to  repeat  the  MoM  computations  in  order  to  ana¬ 
lyze  the  scattering  from  many  different  incidence  angles, 
it  is  desirable  to  decouple  the  computations  of  the  two 
methods.  There  are  two  ways  of  doing  this.  The  first 
method  preserves  the  coupling  interactions  between  dif¬ 
ferent  slots;  the  second  involves  an  approximation  which 
neglects  the  coupling  between  different  slots  to  achieve 
lower  computational  complexity. 

To  decouple  the  MoM  computations  from  the  SBR 
computations  while  preserving  the  coupling  between  the 
various  slots,  the  incident  magnetic  field  on  the  slot  aper¬ 
tures  can  be  expanded  in  terms  of  basis  functions.  As¬ 
suming  that  the  width  of  a  slot  is  much  less  than  its 
length,  the  component  of  the  incident  magnetic  field 
along  the  length  of  a  slot  is  the  only  component  of  in¬ 
terest.  A  convenient  basis  set  is  the  set  of  pulse  basis 
functions,  where  each  function  is  defined  to  be  one  on  a 
portion  of  a  single  slot  aperture  and  zero  elsewhere.  The 
magnetic  currents  on  each  slot  aperture  are  then  com¬ 
puted  with  the  incident  field  on  the  slot  array  set  equal 
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(a)  Geometry 


(b)  Magnitude 


(c)  Phase 


Figure  2:  Comparison  between  the  proposed  MoM/SBR  technique  and  a  previously  published  FEM/SBR  technique. 
The  scattering  results  are  computed  at  a  frequency  of  8  GHz,  and  the  scattering  contribution  of  the  slot  only  is 
shown. 


to  each  of  the  basis  functions  in  turn.  A  matrix- vector 
multiply  is  then  carried  out  during  the  SBR  computa¬ 
tions.  This  matrix- vector  multiply  converts  the  incident 
magnetic  fields  on  the  slot  apertures  to  the  equivalent 
currents  on  the  apertures. 

The  second  method  of  decoupling  the  MoM  computa¬ 
tions  from  the  SBR  computations  neglects  the  coupling 
between  the  individual  slots.  One  slot  on  the  array  is 
chosen,  and  it  is  assumed  that  this  slot  is  the  only  one 
present.  The  MoM  computation  is  carried  out  with  a 
magnetic  field  of  unit  amplitude  on  the  chosen  slot,  and 
the  result  is  a  magnetic  current  on  the  aperture.  It  is 
then  assumed  that  all  of  the  slots  in  the  array  are  equiva¬ 
lent;  the  magnetic  current  on  each  one  is  set  equal  to  the 
incident  magnetic  field  times  the  single  magnetic  current 
computed  by  the  MoM.  This  approximation  significantly 
reduces  the  computational  complexity  and  the  sizes  of 
data  files.  However,  it  does  not  produce  accurate  results 
when  the  frequency  is  near  the  working  frequency  of  the 
slot  array.  This  is  demonstrated  in  Section  4. 


3  Testing 

Before  using  any  new  numerical  technique,  the  technique 
should  be  tested  against  measurements  or  alternate  tech¬ 
niques  to  ensure  its  validity.  Ideally,  measured  data  is 
used  for  comparison,  but  unfortunately  measured  data 
for  this  problem  is  unavailable.  However,  the  validity  of 
the  MoM  computation  involving  the  coupling  between 
the  different  slots  in  the  array  is  validated  by  comparison 
with  previous  MoM  and  finite-element  method  (FEM) 
techniques  [2,7].  The  SBR  method  is  also  validated 
through  extensive,  previous  testing  [4,5].  The  hybrid 
technique  is  validated  by  comparison  with  a  previous 
hybrid  method  to  compute  the  scattering  from  complex 
targets  with  cracks  and  cavities  on  their  surfaces  [3] .  The 
comparison  is  accomplished  by  considering  a  waveguide 
with  each  end  terminated  by  a  (short-circuit)  PEC  plate 
and  with  a  single  slot  on  the  waveguide  surface.  The 
geometry  is  shown  in  Figure  2a.  The  waveguide  with  a 
single  slot  can  be  modeled  both  as  a  slotted  cavity  for 
the  FEM  technique  and  as  a  slotted  waveguide  for  the 
MoM  technique.  Figures  2b  and  2c  show  a  comparison  of 
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Figure  3:  Configuration  of  the  slots  on  the  surface  of  waveguides. 


the  magnitude  and  the  phase  of  the  backscattered  field 
from  the  slot  in  the  presence  of  the  waveguide.  There  is 
good  agreement  both  in  magnitude  and  in  phase.  Both 
magnitude  and  phase  agreement  are  important  so  that 
when  the  scattered  field  from  the  external  geometry  is 
added,  correct  results  are  produced. 

4  Numerical  Examples 

To  show  the  capability  and  utility  of  the  proposed  tech¬ 
nique,  several  numerical  results  are  presented.  For  all 
of  the  numerical  examples,  the  slot  array  contains  16 
waveguides  with  16  slots  on  each  waveguide,  and  the  ar¬ 
ray  is  designed  to  radiate  at  9.1  GHz.  In  addition,  the 
following  parameters  apply  to  all  of  the  examples  pre¬ 
sented:  the  upper  waveguide  wall  in  which  the  slots  are 
cut  is  0.08  cm  thick,  the  waveguides  are  separated  by 
walls  0.1  cm  thick,  each  slot  is  1.6  cm  long  and  0.16 
cm  wide,  and  the  slots  are  positioned  on  the  waveguide 
surface  as  shown  in  Figure  3,  where  the  offset  of  each 
slot  from  the  center  of  the  waveguide  is  0.15  cm.  Unless 
otherwise  noted,  the  coupling  between  individual  slots 
in  the  array  is  included  in  the  results. 

The  first  example  is  a  planar  slot  array  sitting  on  a 
simple  ground  plane.  The  waveguides  are  2.230  cm  wide 
by  1.016  cm  high,  the  slot  centers  are  2.444  cm  apart, 
and  the  first  and  last  slot  centers  are  1.222  cm  from 
the  ends  of  the  waveguides.  Thus,  the  entire  slot  array 
and  the  ground  plate  are  37.3  cm  wide  by  39.1  cm  long. 
In  Figure  4,  the  RCS  of  the  plate  with  the  slots  is  su¬ 
perimposed  on  the  scattering  from  the  plate  alone.  The 
scattering  frequency  is  9.1  GHz,  which  is  the  working  fre¬ 
quency  of  the  slot  array.  Figure  4  shows  results  in  both 
the  H- plane  and  the  .E-plane  and  for  waveguides  which 
are  terminated  both  with  matched  loads  and  with  short 
circuits.  For  matched  waveguide  loads,  each  waveguide 
is  terminated  with  an  impedance  sheet  which  is  matched 
to  the  characteristic  impedance  of  the  guide.  In  the  case 
of  short  circuit  waveguide  loads,  each  waveguide  is  ter- 


mintated  with  a  PEC  plate.  For  some  incidence  angles, 
the  slot  array  has  a  dominant  effect  on  the  scattering. 

The  second  example  is  a  slot  array  on  a  cylinder  with 
a  nose  cone.  The  radius  of  the  cylinder  is  16.096  cm,  and 
the  length  without  the  nose  cone  is  100  cm.  The  nose 
cone  is  30  cm  long.  The  waveguide  cross-sections  are  sec¬ 
toral  in  shape  and  are  1.016  cm  thick.  Along  the  slotted 
surface,  the  waveguides  are  2.230  cm  wide.  The  slots  are 
2.573  cm  apart,  and  the  first  and  last  slots  are  1.287  cm 
from  the  ends  of  the  waveguides.  The  entire  slot  array 
is  37.3  cm  along  the  circumference  of  the  cylinder  and 
41.2  cm  along  the  axis  of  the  cylinder.  In  Figure  5,  the 
if -plane  RCS  of  the  cylinder  alone  and  the  RCS  of  the 
cylinder  with  the  slot  array  are  compared.  The  scatter¬ 
ing  frequency  is  9.1  GHz,  the  working  frequency  of  the 
slot  array,  and  again,  there  are  scattering  directions  for 
which  the  slot  array  dominates  the  return. 

The  next  example  is  intended  to  show  the  effect  of 
the  uncoupled  slot  approximation  which  was  discussed 
in  Section  2.3.  Figure  6  shows  the  RCS  of  the  same 
geometry  considered  in  the  second  example,  but  as  a 
function  of  frequency.  The  incident  direction  is  40°  in 
the  E-plane.  The  RCS  computed  considering  the  cou¬ 
pling  between  individual  slots  is  plotted  with  the  RCS 
computed  by  neglecting  the  slot  coupling.  The  approx¬ 
imation  neglecting  slot  coupling  is  reasonably  accurate 
away  from  the  working  frequency  of  the  slot  array  an¬ 
tenna,  but  there  is  significant  error  near  the  working  fre¬ 
quency.  Thus,  this  approximation  must  be  applied  with 
care. 

The  final  example  shows  the  usefulness  of  the 
method.  The  planar  slot  array  antenna  from  the  first 
example  is  mounted  in  the  nose  cone  of  an  aircraft.  Fig¬ 
ure  7a  shows  the  aircraft,  and  Figures  7b  and  7c  show  the 
HH-polarized  range  profile  of  the  airplane  both  with  and 
without  the  slot  array.  The  range  profile  is  the  time  do¬ 
main  response  to  an  incident  sine  pulse.  The  sine  pulse 
in  this  example  has  a  center  frequency  of  10  GHz  and 
a  bandwidth  of  4  GHz,  and  the  slot  array  has  shorted 
waveguide  loads.  The  slot  scattering  has  an  impact  on 
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(a)  H- plane,  matched  waveguide  loads  (b)  .E-plane,  matched  waveguide  loads 


(c)  E-plane,  short-circuit  waveguide  loads  (d)  E-plane,  short-circuit  waveguide  loads 


Figure  4:  RCS  of  a  planar  slot  array  on  a  ground  plate  at  9.1  GHz,  the  working  frequency  of  the  slot  array. 


the  range  profile. 

5  Conclusion 

A  hybrid  MoM/SBR  method  is  developed  to  compute 
the  scattering  from  a  complex,  3-D  target  with  a  slotted 
waveguide  array  antenna.  Because  the  target  is  large  and 
3-D,  the  MoM  alone  cannot  efficiently  compute  the  scat¬ 
tering,  and  because  the  slots  on  the  waveguides  are  small 
features,  the  SBR  method  alone  is  not  accurate.  The 
hybrid  method  combines  the  two  individual  methods  in 
such  a  manner  that  the  scattering  can  be  efficiently  and 
accurately  computed.  In  the  hybrid  method,  the  MoM 
is  used  to  model  the  details  of  the  slot  array,  and  the 
SBR  method  is  used  to  model  the  electromagnetic  inter¬ 
actions  with  the  large,  complex  target.  The  method  is 
validated  by  comparison  to  previously  published  meth¬ 
ods.  Numerical  examples  show  the  need  to  include  a 
slot  array  model  when  computing  the  scattering  from 
a  complex  target  with  a  slotted  waveguide  array.  The 


examples  also  illustrate  the  capability  of  the  method. 
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49 


-500  -400  -300  -200  -100  0  100  200  300  400  500 

Down  Range  (Inch) 
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(c)  Range  profile  with  slot  array. 
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Effects  of  Gaps  Among  Panels  in 
Radio  Astronomy  Reflector  Antennas 


Giuseppe  Pelosi,  Roberto  Coccioli,  Alessio  Gaggelli 


Abstract —  The  main  reflector  of  antennas  used  for 
radio  astronomy  consists  of  hundreds  of  panels  among 
which,  for  various  reasons,  small  gaps  are  left.  In  this 
paper,  the  effects  of  these  gaps  on  the  field  scattered 
by  the  reflector  are  analyzed  by  means  of  an  hybrid 
numerical  technique  which  combines  the  Finite  Ele¬ 
ment  Method  (FEM)  and  the  Method  of  Moments 
(MoM).  Numerical  results  pertaining  to  the  case  of 
an  incident  plane  wave  are  presented,  and  the  effects 
of  the  introduction  of  corrugations  inside  the  gaps  to 
minimize  the  power  flowing  through  the  gaps  them¬ 
selves  are  discussed. 

I.  Introduction 

THE  construction  of  the  main  reflector  is  the 
most  difficult  mechanical  problem  encountered 
in  the  realization  of  both  radio  telescope  and  large 
antennas  such  as  those  used  for  telecommunications 
in  deep  space.  The  main  reflector  is  realized  as  a  col¬ 
lection  of  conductive  panels  whose  shapes  are  usually 
trapezoidal  or  hexagonal.  The  panels  are  left  sepa¬ 
rated  by  a  small  gap  to  account  for  thermal  expan¬ 
sion,  deformations  due  to  gravity,  and  to  ease  the 
construction  of  the  reflector  itself.  For  radio  astron¬ 
omy  applications,  the  main  issues  are  the  antenna 
maneuverability  (from  a  mechanical  point  of  view) 
and  the  antenna  gain  over  a  wide  frequency  band 
(from  an  electromagnetic  point  of  view). 

The  reflector  for  these  antennas  is  usually  designed 
assuming  its  surface  as  smoothly  curved  and  contin¬ 
uous  and  using  high  frequency  techniques.  This  ap¬ 
proach  is  most  useful,  but  unfortunately  it  is  not  able 
to  account  for  the  presence  of  gaps  whose  effects  in 
the  antenna  performance  are  usually  neglected  [1]. 
More  flexible  and  accurate  numerical  methods,  such 
as  the  Finite  Element  Method  (FEM)  [2]-[3],  can  be 
used  to  simulate  the  presence  of  the  gaps.  Elec¬ 
tromagnetic  scattering  from  gaps  in  perfectly  con¬ 
ducting  planes  has  been  analyzed  by  Gedney  and 
Mittra  [4]  and  Jin  and  Volakis  [5]  using  a  hybrid 
FEM/Method  of  Moments  and  FE/Boundary  Inte- 
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Fig.  1.  Geometrical  model  of  a  gap  between  two  tilted  panels. 

gral  technique,  respectively.  The  aim  of  this  work 
is  to  extend  the  numerical  analysis  to  the  practical 
problem  of  gaps  among  tilted  metallic  planes  as  it  is 
the  case  in  reflector  antennas. 

To  develop  an  efficient  analysis  tool,  the  hybrid 
technique  proposed  in  [6]  to  analyze  electromagnetic 
scattering  from  a  wedge  with  cavity-backed  apertures 
in  its  faces,  is  modified  and  applied  to  this  case.  Re¬ 
sorting  to  this  hybrid  method  limits  the  applications 
of  the  FEM  to  the  gap,  while  using  high  frequency 
techniques  to  simulate  the  large  reflector.  Moreover, 
the  flexibility  of  the  FEM  allows  considering  the  ac¬ 
tual  geometric  configuration  of  the  problem,  in  which 
the  angle  formed  by  two  panels  of  the  main  reflector 
varies  according  to  their  location  in  the  reflector  it¬ 
self,  or  the  panels  present  corrugation  in  their  edges 
to  reduce  transmission  through  the  gap  [7], 

The  hybrid  technique  employed  relies  on  a  particu¬ 
lar  formulation  of  the  equivalence  principle  which  im¬ 
plies  covering  the  gap  apertures  with  metallic  plates. 
The  original  configuration  is  thus  subdivided  into 
two  exterior  problems,  each  of  which  is  comprised  of 
a  perfectly  conducting  wedge,  and  an  interior  prob¬ 
lem,  consisting  of  a  cavity  with  the  same  shape  as 
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the  gap.  The  exterior  canonical  problems  are  for¬ 
mulated  by  exploiting  the  exact  perfectly  conduct¬ 
ing  wedge  Green’s  function,  while  the  more  involved 
interior  problem  is  treated  by  means  of  the  FEM. 
The  solutions  are  then  joined  by  enforcing  the  con¬ 
tinuity  condition  of  the  fields  at  each  gap  aperture. 
This  procedure  restricts  the  application  of  the  FEM 
to  the  small  geometric  domain  constituted  by  the 
gap.  Furthermore,  using  the  Finite  Element  Method 
to  analyze  the  interior  problem  allows  treating  odd 
shaped  gaps,  such  as  those  existing  between  tilted 
panels,  which  cannot  be  easily  taken  into  account 
with  other  numerical  techniques. 

The  hybrid  technique  employed  is  briefly  outlined 
in  Section  II.  The  numerical  solution  implemented 
is  described  in  Section  III.  Finally,  numerical  results 
are  presented  in  Section  IV  to  demonstrate  both  the 
accuracy  of  the  technique  and  the  effects  of  the  gaps. 

II.  Formulation 

Since  the  wavelength  is  always  much  smaller  than 
the  panel  dimensions,  the  reflector  surface  can  be  lo¬ 
cally  approximated  by  a  geometric  canonical  model 
constituted  by  two  half-planes  of  finite  thickness 
forming  an  angle  mr  (Fig.  1).  This  canonical  model 
is  assumed  as  uniform  along  the  z-axis,  which  is  par¬ 
allel  to  the  half-plane  edges,  while  <j>*  and  j3'  denote 
the  direction  of  propagation  of  the  incident  field.  The 
following  analysis  has  been  carried  out  by  consider¬ 
ing  the  case  ft  =  7r/2  and  an  incident  plane  wave 
with  unitary  magnetic  field  parallel  to  the  z-axis 
(TE*  case,  Ez2nc  =  0).  Since  the  evaluation  of  power 
losses  due  to  transmission  of  the  field  through  the 
gap  is  the  parameter  of  interest,  only  this  polariza¬ 
tion  needs  to  be  considered.  In  fact,  since  the  width 
of  the  gap  is  much  smaller  than  the  wavelength  over 
the  frequency  range  of  interest,  TM~  waves  do  not 
propagate  through  the  gap.  The  time  dependence 
exp (jcot)  has  been  assumed  and  suppressed  in  the 
following  analysis. 

The  procedure  proposed  in  this  paper,  which  is 
an  extension  of  that  presented  in  [6],  consists  of 
two  consecutive  steps.  In  the  first  step,  two  ficti¬ 
tious  surfaces  Fi  and  F2  are  defined  by  extending 
in  free  space  the  faces  of  the  perfectly  conducting 
panels,  so  that  they  intersect  at  Q  and  Q/  respec¬ 
tively  (Fig.  2).  Then,  the  equivalence  theorem  is 
applied  so  as  to  subdivide  the  problem  into  three 
separate  configurations  (Fig.  3),  simpler  than  the 
original  one.  Accordingly,  the  surfaces  Ti  and  T2 
are  closed  with  a  perfectly  conducting  cover  and  the 
three  new  problems  are  coupled  through  the  equiva¬ 
lent  surface  magnetic  current  densities  Mext1,2  and 
M*nn’2,  which  are  impressed  on  opposite  sides  of  the 
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Fig.  2.  Surfaces  Ti  and  P2  used  for  the  application  of  the 
equivalence  principle. 

surfaces  Vi#  •  Consequently,  the  original  problem  is 
reduced  to:  z)  the  analysis  of  the  field  in  the  exterior 
region  of  two  complementary  perfectly  conducting 
wedges  (Fig.  3a);  ii )  the  analysis  of  the  field  inside 
the  equivalent  closed  cavity  (Fig.  3b). 

For  case  i),  the  first  exterior  problem  refers  to 
the  front  part  of  the  reflector  and  it  is  necessary 
to  compute  the  magnetic  field  due  to  the  incident 
plane  wave  and  the  unknown  exterior  magnetic  cur¬ 
rent  Mextl  at  Ti  (Fig.  3a),  defined  as  Mextl  =  Exn, 
where  n  is  the  outward  normal  unit  vector  to  T\  and 
E  is  the  electric  field  at  the  same  surface.  The  second 
exterior  problem,  similar  to  the  previous  one  except 
for  the  absence  of  excitation  from  the  backside  of  the 
reflector,  consists  of  the  evaluation  of  the  field  radi¬ 
ated  by  the  equivalent  currents  Mext2  =Exn.  These 
latter  radiate  in  the  presence  of  the  reconstructed 
perfectly  conducting  wedge  with  its  edge  in  Q'  (Fig. 
3a).  Solution  of  the  problem  ii)  implies  the  compu¬ 
tation  of  the  field  radiated  by  the  unknown  magnetic 
currents  Mm*1,2  inside  the  equivalent  closed  cavity 
sketched  in  Fig.  3b. 

The  second  step  of  the  procedure  consists  of 
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b) 

Fig.  3.  a)  Exterior  problems,  b)  Interior  problem.  The 
three  problems  are  coupled  through  the  equivalent  mag¬ 
netic  currents. 

matching  the  exterior  and  interior  problems.  This 
can  be  achieved  by  enforcing  the  continuity,  across 
the  surfaces  r^,  of  the  tangential  component  of  the 
total  electric  and  magnetic  fields.  By  choosing  the 
interior  current  as  Mintl’2  =  -Mext1-2,  the  continu¬ 
ity  of  the  tangential  component  of  the  electric  field  is 
automatically  guaranteed,  while  the  boundary  con¬ 
dition  on  the  magnetic  field  must  be  explicitly  en¬ 
forced: 

n  x  Hextl  =  n  x  H*ntl,  (la) 

n  x  Hext2  =  n  x  Hint2.  (lb) 

Regardless  of  the  matching  procedure  adopted,  the 
exterior  magnetic  field  Hextl  =zHfu  in  equation 
(la),  which  is  defined  in  problem  i ),  may  be  ex¬ 
pressed  as  the  superposition  of  three  contributions: 

Hf1  =  Elnc  +  Hfot  +  Rrzadl(Mextl),  (2) 

where  H*n‘ =  is  the  incident  field,  Wzcat  is  the  field 
scattered  by  the  equivalent  wedge  and  Hxadl(Mexil) 


is  the  field  radiated  by  the  exterior  equivalent  sur¬ 
face  currents  impressed  on  IY  In  the  second  exte¬ 
rior  problem  just  one  term  is  present,  that  is  the 
H ft2  =H:od2(Mext2).  The  scattered  field  Hfat  is 
evaluated  by  using  the  perfectly  conducting  wedge 
Green’s  function  expressed  in  terms  of  eigenfunc¬ 
tions  [8],  or  its  Uniform  Geometrical  Theory  of  Dif¬ 
fraction  (UTD)  approximation  [9],  depending  on 
the  distance  between  the  observation  point  and 
the  edge  Q  of  the  equivalent  wedge.  The  evalu¬ 
ation  of  the  contribution  Hxadl’2(Mext1’2),  due  to 
the  equivalent  magnetic  current  distribution,  has 
been  done  only  with  perfectly  conducting  wedge 
Green’s  function  because  the  surfaces  are  al¬ 
ways  close  to  the  edge  of  the  wedge.  The  interior 
magnetic  field  Hintl>2  =zHfn>2  in  the  right  hand 
side  of  equations  (la)  and  (lb)  is  produced  by  the 
interior  equivalent  surface  current  distribution  im¬ 
pressed  on  H*nt1’2  =H*ntl>2(Mintl,  Mint2)  = 
jjinti,2(  ^jexti  ^  anc[  js  computed  via  the 

FEM  solving  the  scalar  Helmholtz  equation  for  this 
field  component. 

Two  different  matching  techniques,  belonging  to 
the  inward-looking  and  outward-looking  [10]  cate¬ 
gories,  could  be  employed.  Since  the  gaps  among 
panels  have  dimensions  much  smaller  than  the  wave¬ 
length  in  the  frequency  range  of  interest,  the  conti¬ 
nuity  condition  in  equations  (la)  and  (lb)  can  be  di¬ 
rectly  solved  by  a  point-matching  formulation  of  the 
MoM,  as  suggested  in  [11],  The  use  of  this  scheme 
makes  it  necessary  to  explicitly  calculate  the  contri¬ 
bution  arising  from  the  interior  problem  to  the  im¬ 
pedance  matrix  of  the  MoM.  To  this  end,  the  FEM  is 
employed  to  evaluate  the  magnetic  field  on  Ti>2,  due 
to  each  basis  function  used  to  expand  the  interior 
magnetic  current  distribution  Mint1'2  =  —  Mext1,2 
on  Fi,2  itself.  This  procedure  belongs  to  the  inward- 
looking  formulation  class,  and  similarly  to  all  the 
procedures  of  this  kind,  it  is  numerically  efficient.  As 
a  matter  of  fact,  it  first  solves  the  interior  problem, 
which  requires  the  inversion  of  a  sparse,  symmetric, 
and  real  (for  lossless  materials  filling  the  cavity)  ma¬ 
trix,  and  afterwards  it  matches  the  interior  and  the 
exterior  problem.  This  latter  task  requires  the  solu¬ 
tion  of  a  system  of  equations  whose  coefficient  matrix 
is  full,  but  smaller  than  that  related  to  the  interior 
problem. 

This  matching  scheme  also  retains  the  major 
shortcoming  of  inward-looking  formulations:  it  suf¬ 
fers  from  the  possible  presence  of  non-physical  res¬ 
onant  solutions,  which  may  arise  when  solving  the 
closed  interior  problem  in  the  lossless  case.  How¬ 
ever,  due  to  the  dimensions  of  the  interior  problem, 
this  is  not  an  issue  in  this  application. 
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III.  Numerical  Implementation 

The  inward-looking  formulation  for  the  problem 
leads  to  the  standard  MoM  matrix  equation 

_j_  Yextl  -j-  Yint2  -f*  Yex£2]  *  V  —  I  (3) 

where  Yintly 2  and  Yext1,2  account  for  the  contribu¬ 
tions  to  the  magnetic  field  at  T1j2  coming  from  the 
interior  and  the  exterior  problem,  respectively.  Fur¬ 
thermore,  in  equation  (3),  I  is  the  excitation  column 
vector,  due  to  the  incident  plus  scattered  fields  at 
the  aperture  opening  Ti ,  while  V  is  the  column  vec¬ 
tor  of  the  unknown  coefficients  of  the  basis  functions 
used  to  expand  the  equivalent  magnetic  current  dis¬ 
tribution.  Careful  attention  must  be  paid  in  choos¬ 
ing  this  set  of  basis  functions.  As  a  matter  of  fact, 
the  transverse  components  of  the  electric  field  are 
singular  [12],  [13]  at  the  two  edges  of  each  of  the 
thick  panels  delimiting  the  gap  and  this  singularity 
is  apparently  also  present  in  the  equivalent  magnetic 
current.  For  this  reason,  the  four  basis  functions  cen¬ 
tered  at  the  ends  of  the  fictitious  surfaces  and  F2 
have  been  chosen  with  a  singular  behavior  of  the  kind 
—1/(2  ffi),  where  p  is  the  distance  from  the  edge  of 
the  wedge.  Indeed,  this  is  the  kind  of  singularity  that 
the  transverse  component  of  the  electric  field  exhibits 
at  the  proximity  of  a  90°  perfectly  conducting  wedge 
[13].  These  four  basis  functions,  as  well  as  all  the 
other,  span  a  subsection  of  Ti  and  T2  and  are  built 
in  such  a  way  that  they  vanish  at  all  the  other  nodal 
points  on  T i  and  P2  belonging  to  the  same  subsec¬ 
tion.  All  the  other  basis  functions  are  polynomials 
which,  on  the  surfaces  Ti  or  r2,  coincide  with  the 
shape  functions  used  for  the  FEM  discretization  of 
the  interior  problem. 

The  magnetic  field  required  to  calculate  each  en¬ 
try  of  the  matrix  Yintl’2  is  obtained  by  applying  the 
finite  element  method  with  second  order  nodal  ele¬ 
ments  to  the  interior  region.  The  FEM  has  been  cho¬ 
sen  because  it  allows  rather  general  equivalent  cavity 
configurations,  such  as  those  arising  in  case  of  arbi¬ 
trarily  tilted  panels,  corrugations  in  the  thick  side  of 
the  panels,  the  presence  of  dielectric  material,  and  so 
on.  It  is  important  to  remark  that  the  construction 
of  this  (N  x  N )  matrix,  which  needs  the  solution  of 
N  interior  problems  (where  N  is  the  number  of  basis 
functions  used  to  expand  all  the  magnetic  equiva¬ 
lent  currents  on  T\  and  P2),  can  be  done  just  once 
for  each  gap  configuration,  regardless  the  direction 
of  the  incident  wave. 

Yextl>2  in  equation  (3)  results  from  the  term 
H™dl,2(Mexn’2),  and  is  evaluated  by  referring  to  the 
equivalent  problem  of  a  wedge  illuminated  by  a  mag¬ 
netic  line  source.  Also  the  construction  of  this  matrix 
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Fig.  4.  Transmission  coefficient  through  a  0.05  cm  wide  gap 
in  a  0.3  cm  thick  conducting  panel  for  different  angles  of 
incidence. 


can  be  done  just  once  for  a  gap  with  a  specified  aper¬ 
ture,  and  subsequently  used  for  any  direction  of  the 
incident  electromagnetic  field. 

IV.  Numerical  Results 

To  validate  the  code,  previously  published  results 
and  measurements  have  been  considered.  Fig.  4 
shows  the  transmission  coefficient  through  a  slot  with 
width  w  =  0.05  cm  in  a  0.3  cm  thick  conduct¬ 
ing  sheet.  Computations  have  been  performed  for 
three  different  values  of  the  incidence  angle,  and  re¬ 
sults  agree  very  well  with  those  reported  in  [4].  It 
should  be  noted  that  the  transmission  coefficient  has 
been  defined  as  the  ratio  between  the  power  flowing 
through  the  gap  and  that  through  the  same  surface 
located  in  the  free  space  due  to  a  plane  wave  im¬ 
pinging  at  <j>  =  90°,  and  this  accounts  for  the  values 
above  unity  in  Fig.  4. 

A  second  comparison  has  been  made  analyzing  a 
more  involved  configuration,  for  which  measured  re¬ 
sults  are  available  in  the  literature  [14].  It  is  com¬ 
prised  of  a  gap  with  length  3.6  mm  between  two  per¬ 
fectly  conducting  panels  5  mm  thick.  Each  panel 
has  a  corrugation  3  mm  wide  and  7.5  mm  deep.  Fig. 
5  shows,  for  a  such  configuration,  the  transmission 
coefficient  versus  frequency,  where  the  transmission 
coefficient  is  defined  as  the  ratio  of  the  power  flow¬ 
ing  through  the  gap  with  the  corrugation  and  with¬ 
out.  When  the  panels  are  aligned  (n  =  1),  the 
measured  transmission  coefficient  from  reference  [14] 
(dots)  agrees  well  with  the  computed  one  (dashed 
line).  It  can  be  seen  that,  at  a  frequency  close  to  the 
resonant  frequency  of  the  corrugations,  the  transmis¬ 
sion  coefficient  becomes  quite  low,  increasing  the  an¬ 
tenna  performance.  Fig.  5  also  shows  the  transmis- 
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Fig.  5.  Transmission  coefficient  through  a  slit  between  two 
panels  5  mm  thick.  Each  panel  presents  a  corrugation 
7.5  mm  deep  and  3  mm  wide.  </>'  =  75°.  Dashed  line: 
n  —  1)  computed  results;  dots:  n  =  1,  measured  results 
[14];  continuous  line:  n  =  5/6,  computed  results. 

sion  coefficient  versus  frequency  when  the  two  panels 
are  tilted  of  an  angle  equal  to  150°  (n  =  5/6,  contin¬ 
uous  line).  The  length  of  the  gap  Ti  is  kept  constant 
to  3.6  mm,  but  since  the  opposite  faces  of  the  pan¬ 
els  are  moved  farther  from  each  other,  the  resonance 
frequency  of  the  corrugated  gap  decrease,  as  shown 
by  the  computed  results. 

Fig.  6  shows  the  amplitude  of  the  magnetic  field 
longitudinal  component  Hz  inside  the  gap  relative  to 
the  case  n  =  5/6  of  Fig.  5.  The  incident  field  has  a 
frequency  of  10  GHz  and  for  such  value  the  corruga¬ 
tion  depth  is  equal  to  A/4.  Hence,  the  corrugations 
create  an  equivalent  magnetic  wall  along  the  thick 
side  of  the  panels,  which  act  as  a  short  circuit  on 
the  magnetic  field.  As  matter  of  fact,  Fig.  6  shows 
maxima  of  magnetic  field  at  the  end  of  the  corruga¬ 
tions  and  minima  at  their  entrance.  The  low  level 
of  the  magnetic  field  at  the  surface  T2  accounts  for 
the  minimum  in  the  transmission  coefficient  through 
the  aperture.  Although  introducing  the  corrugation 
is  effective  only  in  a  narrow  band,  and  thus  it  is  not 
possible  to  optimize  the  radio  astronomy  reflector 
antenna  performance  in  the  whole  frequency  band  of 
interest,  corrugations  may  be  used  to  optimize  the 
antenna  in  specific  relatively  narrow  bands  as,  for 
instance,  those  also  used  for  telecommunication  pur¬ 
poses. 

The  third  configuration  considered  is  comprised  of 
a  6  mm  wide  gap  between  2  mm  thick  panels  tilted 
at  an  angle  of  170°.  These  parameters  have  been 
measured  directly  on  the  parabolic  radio  astronomy 
reflector  antenna  of  the  Italian  National  Research 
Council  located  at  Medicina  (Bologna,  Italy). 

In  Fig.  7,  the  equivalent  magnetic  currents  at  the 
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Fig.  6.  Amplitude  of  the  magnetic  field  longitudinal  compo¬ 
nent  inside  a  corrugated  gap  between  two  panels  5  mm 
thick  and  tilted  at  150°.  Each  panel  presents  a  corruga¬ 
tion  7.5  mm  deep  and  3  mm  wide.  d/  =  75°  /  =  10 
GHz. 
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Fig.  7.  Normalized  amplitude  of  equivalent  magnetic  current 
distributions  at  surfaces  Fi  and  T2  for  a  slit  6  mm  wide 
and  2  mm  deep  between  two  panels  tilted  at  170°.  The 
incident  plane  wave  impinges  at  an  angle  <pr  =  85°  and 
has  a  frequency  equal  to  10  GHz. 


apertures  Ti  and  T2  are  plotted.  The  incident  plane 
wave  has  a  frequency  equal  to  10  GHz  and  impinges 
at  an  angle  of  <j/  =  85°.  At  both  the  apertures,  the 
equivalent  current  distribution  is  almost  uniform  in 
the  central  part,  while  it  exhibits  the  expected  singu¬ 
larity  near  the  edges  of  the  panels.  These  equivalent 
magnetic  currents  provide  an  almost  isotropic  radi¬ 
ation  pattern.  This  means  that  radiation  incoming 
from  the  backside  of  the  reflector  is  received  through 
the  gap,  thus  increasing  the  antenna  noise  tempera¬ 
ture. 

Further  interesting  results  pertains  to  the  field 
scattered  by  the  two  panels.  Fig.  8  shows  the  nor¬ 
malized  amplitude  of  the  scattered  field  at  a  distance 
of  p  =  10 A  from  the  reconstructed  wedge  Q  versus 
the  observation  angle.  The  two  panels  are  illumi¬ 
nated  by  a  plane  wave  impinging  again  at  <p'  =  85° 
but  with  a  frequency  equal  to  50  GHz.  The  dotted 
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Fig.  8.  Normalized  amplitude  of  the  magnetic  field  at  a  dis¬ 
tance  p  =  10A  from  the  reconstructed  wedge  Q  of  two 
panels  tilted  at  170°.  <f>  =  85°;  /  =  50  GHz.  Dots: 
continuous  panels,  the  presence  of  the  gap  is  neglected. 
Continuous  line:  the  presence  of  a  slit  6  mm  wide  and  2 
mm  deep  is  considered. 


line  refers  to  the  field  computed  by  neglecting  the 
presence  of  the  gap.  This  field,  denoted  as  uo ,  rep¬ 
resents  a  zero-th  order  solution.  The  first  order  ad¬ 
ditive  contribution  u\  is  given  by  the  field  radiated 
by  the  equivalent  magnetic  current  computed  with 
the  technique  proposed.  The  field  uq  +  u\  is  shown 
in  Fig.  8  with  a  continuous  line.  It  is  seen  that  in 
the  backscattering  direction  there  are  no  differences 
between  the  two  patterns,  while  the  contribution  u\ 
of  the  gap  strongly  affects  the  radiation  pattern  in 
the  other  directions,  and  this  effect  increases  with 
frequency. 


V.  Conclusion 

The  effects  of  gaps  among  panels  comprising  the 
main  reflector  of  large  antennas  used  for  radio  as¬ 
tronomy  have  been  studied  by  means  of  a  hybrid  nu¬ 
merical  technique  which  combines  the  FEM  and  the 
MoM.  The  technique  is  based  on  a  particular  appli¬ 
cation  of  the  equivalence  principle  which  allows  the 
reduction  of  the  original  configuration  to  three  sim¬ 
pler  problems:  two  exterior  canonical  problems  and 
one  interior  problem.  The  first  two  can  be  treated  by 
exploiting  the  exact  Green’s  function  or  its  approx¬ 
imation  in  the  framework  of  the  Uniform  Geometri¬ 
cal  Theory  of  Diffraction,  the  latter  can  be  efficiently 
solved  using  the  FEM.  Due  to  the  flexibility  of  the 
FEM,  it  is  possible  to  analyze  arbitrary  configura¬ 
tions  in  an  attempt  to  find  optimized  gap  configura¬ 
tions  able  to  improve  the  antenna  performance. 

The  code  realized  in  this  work  represents  a  first 
step  toward  a  more  sophisticated  tools  able  to  evalu¬ 
ate  the  contribution  to  the  phase  error  at  the  antenna 
opening  due  to  the  presence  of  gaps  between  panels. 


This  contribution,  as  shown  by  the  results  presented, 
may  be  significant. 
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ABSTRACT 

A  numerical  method  based  on  use  of  the  NEC2  code 
has  been  used  to  study  radiation  from  a  slot  in  the 
surface  of  a  conducting  cylinder  and  the  results  used 
in  a  comparison  with  the  eigenfunction  series  solu¬ 
tion.  The  two  cases  of  an  axial  and  circumferential 
slot  have  been  taken  as  representative  between  them 
of  the  more  general  case  of  an  inclined  slot.  In  the 
case  of  a  circumferential  slot,  truncation  of  the  cylin¬ 
der  to  finite  length  has  been  found  in  all  cases  to  lead 
to  numerical  and  analytic  results  which  are  in  poor 
agreement  and  this  has  been  accounted  for  with  an 
argument  based  on  the  geometrical  theory  of  dif¬ 
fraction.  However  this  alone  is  not  enough  to  explain 
what  is  observed  with  axial  slots  where,  for  small  cyl¬ 
inders,  agreement  of  the  two  approaches  is  good  but 
fails  for  larger  cylinders.  It  is  postulated  that  this  is 
due  to  excitation  of  waveguide  modes  on  the  inside  of 
the  larger  cylinders  and  the  result  confirmed  by  clos¬ 
ing  die  ends  of  the  cylinder  in  the  NEC2  model  with  a 
system  of  radial  wires,  when  agreement  is  again  re¬ 
stored. 

1.  INTRODUCTION 

Analytical  studies  of  radiation  from  slots  in  infinite  con¬ 
ducting  cylinders  have  been  made  by  a  number  of  au- 
thors[l],  [2],  [3],  [4],  [5].  However,  when  as  must  neces¬ 
sarily  be  the  case  in  practice,  the  cylinder  in  which  the 
slot  is  cut  is  of  finite  length,  to  find  the  radiation  pattern, 
resort  has  to  be  made  to  either  numerical  methods  or 
some  approximate  technique  of  which,  for  a  sufficiently 
large  cylinder,  the  geometric  theory  of  diffraction  might 
be  an  example. 

This  paper  presents  some  experiences  with  a  numerical 
means  of  determining  the  radiation  pattern  of  a  narrow 
slot  cut  in  the  (assumed  infinitesimally  thin)  wall  of  a 
hollow,  perfectly  conducting,  finite,  circular  cylinder. 
The  application  which  we  have  in  mind  is  personal  sat¬ 
ellite  wireless  mobile  communications  and  our  interest 
will  be  with  narrow  slots,  the  length  of  which  are  no 
more  than  comparable  with  the  half  wavelength  (reso¬ 
nant  or  near  resonant  slots). 


The  paper  is  arranged  as  follows.  Section  2  reproduces 
the  eigenfunction  expansion  of  the  electric  field  ra¬ 
diated  by  an  arbitrarily  oriented  slot  cut  in  an  infinite 
cylinder.  Next  the  radiation  pattern  of  a  slot  cut  in  an 
otherwise  similar  finite  cylinder  is  generated  using  the 
NEC2  package  and  compared  with  the  eigenfunction 
solution  and  in  the  following  section  there  is  a  discus¬ 
sion  of  the  notable  points  of  difference  and  the  reasons 
for  them.  The  paper  concludes  with  a  summary  of  its 
findings. 

2.  AN  ARBITRARILY  ORIENTED  SLOT  ON  AN 
INFINITE  CYLINDER 

When  the  cylinder  is  of  infinite  length,  an  analytical 
solution  for  radiation  from  an  arbitrarily  oriented  slot  is 
possible  in  terms  of  cylindrical  wave  functions  and  has 
been  known  for  a  considerable  time  [4].  Figure  1  shows 
a  coordinate  system  in  terms  of  which  it  is  convenient  to 
write  down  the  eigenfunction  series  expansion.  The 
centre  of  the  slot  lies  in  the  XOZ  coordinate  plane  with 
OX  passing  through  its  centre,  which  is  at  (a, 0,0),  a  be¬ 
ing  the  radius  of  the  cylinder.  The  slot  is  excited  by  a 
voltage  V  applied  across  it. 


Figure  1  Geometrical  parameters 

With  respect  to  local  coordinates  (r|,  Q  drawn  along  the 
centrelines  of  the  slot,  it  is  assumed  that  its  electric  field 
distribution  is  oriented  parallel  to  the  T|  direction,  in 
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which  it  is  also  uniform,  while  in  £  it  is  sinusoidal.  For 
this  configuration,  the  series  solutions  for  the  two  spher¬ 
ical  components  of  the  radiation  field  are 

-  Vocosxpe-j*  v  in+le-^InPn 
te  ~  2n2r sin0  n^B2\kasmd)  1  ’ 
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\|/  is  the  inclination  angle,  l  is  the  slot  length,  w  is  the 
slot  width  and  k  is  the  free  space  wavenumber. 


Getting  to  the  essentials  of  the  problem  is  best  facilitated 
by  considering  the  two  special  cases  \\f  =  0  and  \jr  =  90 , 
when  the  slot  is  respectively  circumferential  and  axial. 
Calculations  for  the  radiation  pattern  in  the  0  =  90  plane 
(the  XOY  plane)  for  these  two  cases  are  shown  in  Fig¬ 
ures  2,  3  and  4. 


(a) 


(b) 


Figure  2.  Comparison  of  normalised  far-field  pattern  in  the 
0=90°  plane  obtained  using  an  analytic  solution  and  NEC2 
for/=2GHz,  fca=0.8,  w=0.05A,,  L =X  and  (a)  an  axial  with 
1= 0.5  A,  and  (b)  a  circumferential  slot  with  0.3X. 


Figure  3.  Comparison  of  normalised  far-field  pattern  in  the 
0=90°  plane  obtained  using  an  analytic  solution  and  NEC2 
for/=2GHz,  1.5,  w=0.05?i,  L=X  and  (a)  an  axial  with 
1=0. 5X  and  (b)  a  circumferential  slot  with  1=0  3X. 


Figure  4.  Comparison  of  normalised  far-field  pattern  in  the 
0=90°  plane  obtained  using  an  analytic  solution  and  NEC2 
for/=2GHz,  te=2,  w=0.05X,  h=X  and  l=0.5X  for  (a)  an 
axial  and  (b)  a  circumferential  slot  respectively. 


3.  NUMERICAL  SOLUTION  OF  THE  FINITE 
CYLINDER 

The  NEC2  package  has  been  used  to  solve  the  same 
problems  for  a  finite  cylinder.  NEC2  uses  a  wire  grid 
approximation  for  the  cylinder  which,  when  cut  and  un¬ 
folded  into  a  developed  view,  is  the  simple  rectangular 
mesh  shown  in  Figure  5.  In  modelling  the  axial  slot,  the 
spacing  between  the  vertical  wires  is  chosen  to  be  the 
desired  slot  width.  The  same  is  done  for  the  horizontal 
wires  in  the  case  of  the  circumferential  slot. 
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Figure  5.  NEC2  Wire  Grid  Model 


Three  different  cylinder  sizes  were  considered  for  each 
slot  orientation.  These  were  ka  =  0.8, 1.5  and  2.0.  In  each 
case  the  cylinder  length  was  chosen  as  one  wavelength. 
Figures  2, 3  and  4  superimpose  the  normalised  (for  each 
pattern,  to  its  own  maximum)  NEC-derived  and  analyt¬ 
ic  solutions.  The  choice  to  compare  normalised  patterns 
was  made  as  it  renders  differences  between  the  two  most 
easily  apparent.  However,  it  is  to  be  remarked  that 
where  there  is  strong  agreement,  as  in  Figures  2a  and  3a 
for  example,  at  maximum,  there  is  less  than  0.5  dB  dis¬ 
placement  between  the  two  absolute  patterns.  Not  sur¬ 
prisingly,  where  the  correlation  is  poor,  as  in  Figures  2b 
and  3b,  differences  between  maxima  of  up  to  3  dB  are  in 
evidence.  It  is  interesting  to  compare  the  normalised 
patterns  to  seek  reasons  for  observed  discrepancies. 


as  radiation  in  XOZ  is  concerned,  looking  at  the  devel¬ 
oped  view  of  the  cylinder  shown  in  Figure  6,  we  see  that 
creeping  waves  which  reach  the  shadow  boundary  by 
making  the  minimum  transit  of  three  quartos  of  a  turn 
around  the  cylinder  -  and  these,  rather  than  the  more  fre¬ 
quently  encircling  creeping  waves  which  will  also  be 
present,  are  those  which  make  the  major  contribution  to 
the  radiation  field  -  can  only  be  present  in  90-a  <  0  < 


90+a  ,  where  a  =  atan 


the  neighbourhood  of 


the  XOY  plane.  More  of  these  creeping  waves  will  be 
included  by  using  a  longer  cylinder. 


2iza 

notes:  SB  refers  to  shadow  boundary 

SPP1  and  SPP2  are  the  stationary  phase  points 

Figure  6.  Developed  view  of  the  slotted  cylinder 


4.  DISCUSSION 

The  geometric  theory  of  diffraction  [6]  provides  a  good 
paradigm  in  which  to  look  for  the  answers.  Especially 
for  the  ka  =  0.8  case,  the  cylinders  which  we  have  cho¬ 
sen  would  lie  at  the  lower  end  of  the  range  for  which  this 
theory  can  be  expected  to  give  accurate  answers  but,  at 
least  for  the  purposes  of  qualitative  reasoning,  that  is  not 
of  great  account. 

In  addition  to  direct  radiation  from  the  slot,  on  the  sur¬ 
face  of  the  cylinder  it  will  also  give  rise  to  a  family  of 
creeping  waves  which  will  travel  over  helical  paths 
(straight  lines  on  a  developed  view),  shedding  tangen¬ 
tial  rays  as  they  progress.  Thus,  for  example,  contribu¬ 
tions  to  the  radiation  in  the  XOZ  plane  will  emanate 
from  the  shadow  boundaries  formed  by  the  lines  of  in¬ 
tersection  of  the  plane  YOZ  with  the  cylinder. 

When  the  cylinder  is  of  finite  length,  these  helical  paths 
will  eventually  intersect  with  its  open  ends,  when  two 
things  will  occur. 

One  is  that  creeping  rays  which  would  have  reached  the 
far  field  point  by  being  shed  from  parts  of  the  cylinder 
now  removed  by  truncation  can  no  longer  do  so.  So  far 


The  other  effect  is  that  there  will  be  an  edge  diffraction 
contribution  from  the  open  ends  of  the  cylinder,  the 
dominant  part  of  which  will  come  from  certain  station¬ 
ary  phase  points  which  can  be  located  by  use  of  Fermat’s 
principle  [6].  For  radiation  into  the  positive  XOZ  half 
plane  (<(>  =  0°)  these  will  be  the  points  of  intersection  of 
that  half  plane  with  the  ends  of  the  cylinder,  immediate¬ 
ly  above  and  below  the  slot.  For  other  values  of  <j>  in  0  = 
90°,  the  stationary  phase  points  will  move  circumferen¬ 
tially  around  the  top  and  bottom  edges  of  the  cylinder  in 
sympathy  with  $  (but  through  an  angle  which  depends 
on  the  length  of  the  cylinder  and  which  is  always  less 
than  (|>). 

The  degree  to  which  these  perturbations  can  be  expected 
to  bring  about  significant  deviation  from  the  infinite  cyl¬ 
inder  result  will  vary  markedly  in  the  case  of  an  axial  as 
against  a  circumferential  slot.  The  reason  for  this  is  not 
hard  to  see.  In  both  cases  the  slot  will  behave  like  a  mag¬ 
netic  dipole  lain  on  the  surface  of  the  cylinder  and  the 
source  pattern  of  a  dipole  is  maximum  transverse  to  its 
axis  and  zero  along  its  length.  In  the  case  of  an  axial  slot, 
this  means  that  not  only  are  the  stationary  phase  points 
weakly  illuminated,  the  more  so  for  longer  cylinders 
and  indeed  not  at  all  for  <|>  =  0°,  but  the  part  of  the  family 
of  creeping  waves  which  is  truncated  off,  being 
launched  at  angles  closer  to  the  dipole  axis,  is  also  weak- 
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er.  With  the  circumferential  slot,  exactly  the  opposite  is 
true. 

Figure  7  shows  the  effect  of  increasing  cylinder  length 
for  the  case  of  ka  =  2.0.  In  the  light  of  the  above  it  is  no 
mystery  that  there  is  poor  agreement  between  the  eigen¬ 
function  series  and  the  NEC2  solution  in  the  case  of  a 
circumferential  slot. 


Figure  7.  Comparison  of  normalised  far-held  pattern  in  the 
0=90°  plane  obtained  using  an  analytic  solution  and  NEC2 
for/=2GHz,  £a=2.  vv=0.05^  /=0.5X  and  L=lX}  2X  and  3X 
for  (a)  an  axial  and  (b)  a  circumferential  slot. 

The  problem,  however,  is  still  not  entirely  solved  for  the 
axial  slot.  It  is  only  with  the  smallest  diameter  cylinders 
that  good  agreement  is  obtained.  As  the  cylinder  be¬ 
comes  larger,  there  must  be  an  expectation  that  both  the 
slot  directly  and  edge  diffraction  at  the  ends  will  excite 
propagating  waveguide  modes  inside  it.  Power  coupled 
into  these  modes  will  be  conveyed  without  attenuation 
to  the  open  ends  of  the  cylinder,  from  where  it  will  be 
radiated  to  interfere  with  the  fields  produced  by  the  vari¬ 
ous  mechanisms  already  identified.  When  this  has  the 
possibility  of  happening  can  be  determined  by  looking 
at  the  cutoffs  of  the  modes  in  cylindrical  waveguide  [7], 

TheTEll  mode  has  the  lowest  cutoff  at  )^a  -  1.841  and 
there  is  no  other  mode,  TE  or  TM,  which  has  a  cutoff  be¬ 
low  kcd  =  2.405.  In  the  case  of  thefaz  =  0.8  cylinder,  we 
are  therefore  well  away  from  any  possibility  of  a  wave¬ 
guide  mode  being  present  to  upset  the  agreement  be¬ 
tween  the  eigenfunction  and  NEC2  solutions.  The  oppo¬ 
site,  of  course,  is  true  with  ka  =  2.0  and  is  held  to  explain 
the  disagreement.  The  ka  =  1.5  case  is  more  interesting. 
Taken  at  face  value,  the  evidence  is  that  the  cylinder  is 
below  waveguide  cutoff,  when  there  should  be  no  prob¬ 
lem.  However,  it  is  close  enough  to  cutoff  that,  with  a 
relatively  short  cylinder,  evanescent  fields  will  still 
reach  the  open  ends  of  the  cylinder,  probably  in  suffi¬ 
cient  strength  to  account  for  some  of  the  minor  discrep¬ 
ancies  seen  in  Figure  3a. 


A  way  of  testing  these  conclusions  is  to  repeat  the  NEC2 
solution  with  a  closed  ended  cylinder,  modelled  in  this 
case  by  a  system  of  radial  wires  across  each  end.  This 
will  leave  nearly  everything  else  in  the  problem  intact 
(save  that  diffraction  at  the  ends  is  now  ff om  a  right 
angle  wedge  rather  than  a  half  plane)  while  eliminating 
any  contribution  from  waveguide  modes.  Figure  8a 
makes  the  point;  directly  any  possibility  of  a  waveguide 
contribution  is  removed,  agreement  in  the  axial  case  be¬ 
tween  the  eigenfunction  and  NEC2  solutions  is  much 
improved. 

None  of  this,  of  course,  is  to  say  that,  in  the  case  of  the 
circumferential  slot,  waveguide  modes  will  not  exist  or 
also  have  the  possibility  of  contributing  to  discrepancies 
between  the  eigenfunction  and  numerical  solutions,  but 
it  is  to  say  that  in  this  instance,  complete  salvation  is  not 
to  be  had  simply  by  their  removal.  As  Figure  8b  shows, 
while  the  situation  is  ameliorated,  large  discrepancies 
remain. 


Figure  8.  Comparison  of  normalised  far-field  pattern  in  the 
0=90°  plane  obtained  using  an  analytic  solution  and  NEC 2 
for/=2GHz,  kcv=l,  w=0.05X,  l=0.5X  and  L=  IX  for  (a)  an 
axial  and  (b)  a  circumferential  slot. 


5.  CONCLUSION 

A  wire  grid  model  based  on  the  NEC2  code  has  been  de¬ 
veloped  with  which  to  simulate  radiation  from  near  res¬ 
onant  slots  cut  in  the  surface  of  a  finite  length  conduct¬ 
ing  cylinder  and  the  results  compared  with  those 
obtained  from  an  associated  analytic  solution  which  as¬ 
sumes  an  infinitely  long  cylinder.  Discrepancies  be¬ 
tween  the  two  have  been  explained  using  an  argument 
based  on  the  geometrical  theory  of  diffraction  which 
lends  credibility  to  the  numerical  solution  and  confirms 
its  usefulness  for  the  study  of  radiation  from  arbitrarily 
oriented  slots,  which  are  believed  to  have  potential  ap¬ 
plication  in  the  area  of  personal  mobile  satellite  wireless 
communications. 
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An  interesting  intervention  in  the  final  result  from  wa¬ 
veguide  modes  propagating  within  the  cylinder  but  not 
intended  to  be  part  of  the  model  has  also  been  observed. 
This  serves  to  emphasise  a  need  constantly  to  be  on 
guard  to  ensure  that  a  model  is  chosen  which  simulates 
only  the  problem  of  interest  and  nothing  more.  This  ob¬ 
servation  may  also  make  some  contribution  to  comment 
which  has  appeared  in  NEC  user-group  e-mail  discus¬ 
sions  (NEC-LIST)  on  the  apparent  inability  of  NEC2 
based  solutions  to  match  those  obtained  analytically  for 
larger  diameter  cylinders. 
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Abstract 

We  consider  a  three-dimensional  approach  based  on  the 
KirchhofTs  method  in  order  to  predict  electromagnetic 
waves  propagation  in  urban  environments.  In  particular, 
we  are  interested  here  in  the  evaluation  of  the  electro¬ 
magnetic  field  on  very  large  three-dimensional  domains 
(typically  with  linear  dimensions  of  the  order  of  hun¬ 
dreds  of  meters)  generated  by  a  high  frequency  source 
(typically  of  the  order  of  1GHz  which  corresponds  to  a 
wavelength  of  about  30  centimeters).  Some  numerical 
tests  and  comparisons  with  experimental  measurements 
have  been  done  to  validate  this  approach. 

1  Introduction 

Efficient  planning  of  mobile  communication  systems  is 
based  on  an  accurate  determination  of  the  coverage  re¬ 
gion  of  the  antenna.  In  urban  environments,  this  can 
be  quite  complicated,  since  electromagnetic  waves  are 
absorbed,  reflected,  transmitted  and  diffracted  by  build¬ 
ings.  The  subject  of  our  research  activity  is  the  analysis 
and  development  of  simulation  tools  to  effectively  predict 
the  propagation  of  the  electromagnetic  waves  and  the 
antenna  coverage  in  three-dimensional  urban  domains 
which  have  huge  dimensions  with  respect  to  the  wave¬ 
length. 

Techniques  used  to  solve  this  problem  in  a  reason¬ 
ably  fast  way  include  ray  tracing  (Ikegami  et  al ,  1991, 
Rossi  et  al ,  1992,  Kiirner  et  al ,  1993,  Rizk  et  al ,  1994) 
and  two-dimensional  Transmission  Line  Matrix  (TLM) 
(Luthi  et  al ,  1996).  We  propose  a  three-dimensional  ap¬ 
proach  based  on  the  Kirchhoff  integral  method.  The 
presence  of  multiple  reflections  has  led  to  implement 
Kirchhoff Js  method  using  an  iterative  procedure  where 
each  iteration  corresponds  to  a  reflection  on  a  surface. 
In  addition,  a  number  of  approximations  have  been  in¬ 
troduced  to  limit  the  computational  complexity. 


This  approach  is  validated  through  numerical  com¬ 
parisons  with  experimental  measurements.  Further  re¬ 
sults  related  to  the  convergence  and  time  scaling  of  the 
method  as  well  as  its  application  to  determine  the  an¬ 
tenna  coverage  in  urban  environments  are  presented  and 
discussed. 

2  The  vector  Kirchhoff  integral 
relation 

In  this  section  we  review  the  theoretical  basis  of  the  vec¬ 
tor  Kirchhoff  integral  relation  (Jackson,  1975,  pages  432- 
435). 

Suppose  that  the  vector  field  E  has  harmonic  time 
dependence  e~iwt .  When  the  field  sources  are  outside 
the  volume  V ,  the  field  E  satisfies  the  vector  Helmholtz 
wave  equation  inside  V , 

(V2  +  fe2)E(x)  =  0  (1) 

where  k  =  ^  and  c  is  the  wave  propagation  velocity.  By 
applying  the  divergence  theorem 

VdV  =  da  (2) 

V  Sv 

it  is  easy  to  see  that,  when  equation  (1)  is  verified,  the 
following  relation  is  an  identity: 

/Akr  Akr 

[=(»'  v's;)-b<n' <3) 

Sv 

In  (3)  r  =  x  -  x' ,  the  point  x  is  inside  the  volume  V 
bounded  by  the  surface  Sv  and  the  unit  normal  n'  is 
directed  into  the  volume  V . 

Equation  (3),  which  is  known  as  the  Kirchhoff  inte¬ 
gral,  expresses  the  electric  field  E  at  a  point  x  inside  a 
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closed  volume  V  in  terms  of  the  values  of  the  field  and 
its  normal  derivative  on  the  boundary  surface  Sv  -  By 
using  Maxwell  equations,  we  can  rewrite  equation  (3)  in 
terms  of  the  electric  and  magnetic  fields  E  and  B: 


E(x) 


j_  /V*r 

47 t  J  r 
Sv 


k  x  (n'  x  E) 


(4) 


+k(ri  x  B)  -  k(n'  •  E) 


dx' 


approximately  into  an  illuminated  region  Si  and  a  shad¬ 
owed  one  Ss .  In  the  shadowed  region  the  total  field  on 
the  surface  Ss  is  nearly  zero  (the  object  is  opaque).  We 
assume  that  there  the  scattered  field  on  Ss  is  equal  and 
opposite  to  the  incident  field.  In  addition,  we  suppose 
that  the  building’s  external  walls  are  either  flat  or  with  a 
radius  of  curvature  large  with  respect  to  the  wavelength. 
As  a  result,  in  the  illuminated  region,  Fresnel  coefficients 
for  reflection  are  used  to  evaluate  the  scattered  field  on 
Si .  This  procedure  can  be  summarized  as  follows: 


For  the  Helmholtz  equation  (1)  to  be  satisfied  inside 
V ,  the  medium  contained  in  V  must  be  homogeneous. 
Therefore,  when  we  need  to  evaluate  the  electromagnetic 
scattering  from  objects  which  have  a  dielectric  constant 
€  i=-  *0j  we  can  identify  V  with  the  free  space  between 
the  scatterers.  The  boundary  surface  Sv  of  V  is  taken 
equal  to  S  U  Sc 0  where  5  is  the  surface  of  the  scattering 
objects  and  Soo  is  the  surface  “at  infinity”.  It  can  be 
proven  (Jackson,  1975,  pages  433-434)  that  the  integral 
over  Soo  vanishes;  so  the  integral  in  (4)  can  be  computed 
only  over  5. 

It  is  useful  to  specialize  (4)  to  a  scattering  situation 
and  to  exhibit  a  formal  expression  for  the  scattering  am¬ 
plitude  as  an  integral  of  the  scattered  fields  over  S.  Since 
the  scattering  objects  are  supposed  to  be  outside  the 
volume  V,  equation  (4)  holds  for  the  scattered  fields 
(ES,BS),  that  is,  the  total  fields  (E,  B)  minus  the  in¬ 
cident  fields  (Ei,Bi)  (i.e.,  the  fields  generated  by  ac¬ 
tive  sources  such  as  antennas).  In  equation  (4)  E(x) 
depends  explicitly  on  the  outgoing  direction  of  k.  The 
dependence  on  the  incident  direction  is  implicit  in  the 
scattered  fields  E  and  B. 

3  Description  of  the  model 


sl)  the  incident  fields  E*  and  B2  (generated  by  the  an¬ 
tenna)  are  evaluated  on  the  surface  S  of  the  scat¬ 
tered 

s2)  the  scattered  fields  E5  and  Bs  on  5  are  approxi¬ 
mated  as  described  above; 

s3)  the  scattered  fields  throughout  the  space  V  are  com¬ 
puted  using  the  Kirchhoff  integral  method. 

This  method  gives  a  good  estimate  of  the  scattered 
field  around  an  opaque  object.  If  the  diffraction  effects 
on  the  building  wedges  are  negligible,  we  can  simplify 
the  above  procedure  by  assuming  a  “Geometrical  Optic¬ 
s’’  shadow  for  the  incident  field,  i.e.  by  assuming  that 
both  the  incident  and  scattered  fields  are  equal  to  zero 
beside  the  obstacle.  When  we  are  dealing  with  a  large 
number  of  scatterers,  this  approximation  can  allow  a  sen¬ 
sible  reduction  in  the  computational  time,  because  it  re¬ 
duces  the  number  of  pair-interactions  (only  the  pairs  of 
scatterers  that  “see”  each  other  interact).  This  aspect 
will  be  discussed  in  detail  in  section  4.1.  Later  on  this 
paper,  we  will  refer  to  the  first  procedure  as  the  “general 
procedure”  and  to  the  second  as  the  “GOS  (Geometrical 
Optics  Shadow)  procedure”. 


The  evaluation  of  the  field  E  through  equation  (4),  re¬ 
quires  the  knowledge  of  the  fields  Es  and  Bs  on  the  sur¬ 
face  5. 

3.1  A  single  scatterer 

In  the  absence  of  knowledge  about  the  correct  fields  E5 
and  Bs  on  the  surface  5,  we  must  make  some  approxi¬ 
mations. 

First,  we  neglect  the  field  transmitted  through  the 
building.  This  is  equivalent  to  assuming  that  the  build¬ 
ing  is  an  “opaque”  object.  This  assumption  is  reasonable 
at  the  frequencies  we  are  dealing  with  as  the  part  of  the 
wave  energy  which  is  not  reflected  from  the  surface,  is 
largely  dissipated  inside  the  building.  Because  the  wave¬ 
length  is  small  with  respect  to  the  linear  dimensions  of 
the  obstacle,  the  surface  S  of  the  building  can  be  divided 


3.2  More  than  one  scatterer 

When  dealing  with  more  than  one  building,  the  proce¬ 
dure  for  one  scatterer  can  be  extended  using  an  iterative 
technique  where  each  iteration  corresponds  to  a  reflec¬ 
tion  on  a  surface. 

a:  The  incident  field  produced  by  all  the  sources  is  eval¬ 
uated  on  each  building  surface  Sn : 

•  at  the  first  iteration,  the  only  source  is  the  an¬ 
tenna; 

•  starting  from  the  second  iteration,  the  sources 
of  incident  field  are  the  antenna  and  all  the  sur¬ 
faces  which  have  been  illuminated  during  the 
previous  iteration.  The  second  contribution  is 
evaluated  by  using  the  Kirchhoff  integral. 
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— —  direct  path 

-  path  involving  one  reflection 

-  path  involving  two  reflections 


Figure  1:  Example  of  source-receiver  paths  involving  up 
to  two  reflections.  These  paths  are  accounted  for  by  the 
method  after  two  iterations 

b:  When  the  incident  fields  E;  and  B ;  are  known  at 
each  point  of  the  surface  5n,  the  scattered  fields  Es 
and  B5  on  Sn  are  evaluated  as  described  in  section 
3.1.  At  the  frequencies  considered,  the  diffraction 
fields  are  small  with  respect  to  the  Geometrical  Op¬ 
tics  fields  and  they  can  be  neglected  if  the  intensity 
of  incident  radiation  is  low.  Therefore,  we  have  cho¬ 
sen  to  evaluate  the  scattered  fields  on  Sn  by  using 
the  “general  procedure”  when  the  source  of  the  inci¬ 
dent  fields  E;  and  B *  is  the  antenna  and  the  “GOS 
procedure”  elsewhere.  E5  and  Bs  will  be  used  in 
the  first  step  of  the  next  iteration  to  compute  the 
incident  field  on  the  other  surfaces. 

c:  When  the  electromagnetic  fields  do  not  significantly 
change  from  one  iteration  to  the  next  one,  the  pro¬ 
cedure  ends. 

d:  At  the  end  of  the  iterative  procedure,  we  know  the 
electromagnetic  fields  on  all  the  involved  surfaces 
and  we  can  use,  once  again,  the  Kirchhoff  integral 
to  evaluate  the  electromagnetic  field  throughout  the 
domain  V . 

The  above  procedure  is  sketched  in  figure  1  in  the  case  of 
two  buildings  and  up  to  the  second  iteration.  It  is  seen 
that,  each  iteration  is  equivalent  to  include  new  paths 
connecting  the  source  to  the  receiver. 

To  analyze  the  convergence  of  the  method,  we  consider 
the  configuration  described  in  figure  2.  A  plane  wave 
of  wavelength  A  =  6.28  meters  and  wave  vector  k  = 
(0.89, 0.42,  —0.15),  whose  projection  onto  z  =  0  is  drawn 
in  the  picture,  has  been  taken  as  the  incident  field.  At 
each  iteration,  the  field  values  are  measured  at  points 


located  at  a  height  of  1.8  meters,  along  the  axis  shown 
in  the  figure. 


Figure  2:  Geometry  used  to  test  convergence  of  the 
method.  Receivers  are  located  along  the  dashed  line  be¬ 
tween  the  two  buildings  (all  dimensions  are  given  in  meters) 

At  iteration  n,  the  convergence  of  the  method  is 
checked  through  the  following  parameters: 

er|;  =  max  {\E?(k)  -  E?~\k)\,k  =  l,Nr}  , 

i  =  1 , 3 

er£  =  max  {er£ . ,  2  =  1,3} 

where  Nr  is  the  number  of  receivers. 

Figure  3  reports  the  values  of  er£  at  each  iteration. 
The  lines  refer  to  two  different  simulations  where  the 
buildings  are  considered  either  metallic  or  built  with  ma¬ 
terial  of  dielectric  constant  e  =  4.  In  the  second  case, 
the  convergence  is  faster  because,  at  each  iteration,  part 
of  the  incident  energy  is  absorbed  by  the  external  walls. 

In  figure  4  we  report  the  values  of  two  components  of 
the  electromagnetic  field,  computed  in  the  case  of  metal¬ 
lic  buildings.  The  convergence  of  the  component  Ey  is 
much  slower  than  that  of  the  component  Ex .  This  can  be 
intuitively  explained  by  noting  that  an  y-  directed  elec¬ 
tric  field  can  be  associated  to  an  x-directed  wave,  which 
is  reflected  between  the  metallic  walls  as  in  a  resonat¬ 
ing  cavity,  while  the  x-directed  electric  field,  related  to 
an  y-directed  wave,  is  free  to  leave  the  system.  So  we 
expect  that  the  y  component  will  keep  changing  at  each 
iteration  (and  further  reflection)  between  the  buildings 
(i.e.  65  <  y  <  100)  more  than  outside  them,  and  this 
is  confirmed  by  the  results  of  the  simulation  shown  in 
figure  4  (bottom). 
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Figure  3:  Convergence  of  the  method.  White  dots  refer  to 
perfectly  conducting  surfaces.  Black  dots  refer  to  dielectric 
surfaces  (  c  =  4  ).  The  quantity  "maximum  variation", 
reported  on  y  axis,  is  defined  as  er^  =  max  {er^.,  i  = 
1,3},  where  er£.  is  er£(  =  max{|£f(fc)  -  E^~1(k)\,  k  = 
1 ,  iVr  }  ,  2  =  1,3  and  Nr  is  the  number  of  receivers 


4  Computational  aspects 

In  this  section  we  present  the  main  computational  fea¬ 
tures  of  the  model. 


4.1  Scaling  of  the  computer  time  with 
the  dimension  of  the  domain 


The  model  described  in  section  3  requires,  on  each  sur¬ 
face,  the  evaluation  of  the  field  scattered  from  the  other 
surfaces.  Thus,  at  first  sight,  the  computational  time 
grows  as  N 2  with  N  the  number  of  surface  elements  and 
as  D 4  with  D  the  diameter  of  the  domain.  This  behavior 
is  even  worse  than  the  one  shown  by  the  methods  which 
require  the  discretization  of  the  whole  space  (D3),  but, 
as  we  will  show,  can  be  easily  improved. 

Consider  the  situation  described  in  figure  5.  S2  is  a 
region  of  the  surface  on  which  we  need  to  compute  the 
field  scattered  by  a  region  S\  of  another  surface;  we  can 
write 

r  =  r0  -  x'  +  x" 


and 

r  =  [rl  -  2r„  •  (x'  -  x")  +  |x'  -  x'f]  * 

When  r  is  large  compared  with  x '  and  x" ,  we  can  rewrite 
r  using  a  suitable  Taylor  expansion: 


r  —  r0  - 


r0  •  (x7  —  x7/)  ]  T 
ro 


where 

|x'-x»[2  [r0  •  (x'  —  x")]2 

2r0  2 rg 


(5) 

(6) 


Ex 


Figure  4:  Components  Ex  (top)  and  Ey  (bottom)  at  dif¬ 
ferent  iterations  in  the  example  shown  in  figure  2.  The 
spatial  coordinate  along  the  dashed  line  is  reported  on  the 
x  axis 


Provided  that  T  is  negligible,  we  have: 


gifcr  _  eikr 0  e-ik*(x/~x//) 


(7) 


Let  d'  and  d"  be  the  diameters  of  the  regions  S\  and 
S2.  To  get  a  quantitative  estimate,  we  impose  that  the 
error  on  the  phase  is  bounded  by  Then,  looking  at 
equation  (6),  we  can  see  that  T  is  negligible  when 


(d‘  +  d")2  A 
8r0  <  16 


(d'  +  d")  < 


(8) 


which  is  equivalent  to  the  well-known  far  field  relation. 
In  this  approximation,  equation  (4)  becomes: 


Es(x) 


i  eikr° 
4? r  ro 


eikx"F(k) 


(9) 
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S2 

Figure  5:  Scattering  geometry 


F(k)  =  ye~ikX  (Kn'  x  B*)  +  tk  x  (n/  x  E") 

Si 

— k(n'  •  E,)]  (l  +  £:)  }  dx>  (10) 

In  equation  (10),  the  integral  does  not  depend  on  the 
coordinate  x"  internal  to  the  region  S2-  Therefore,  when 
we  decompose  the  surfaces  into  regions  which  are  small 
compared  to  their  relative  distance,  we  need  to  evaluate 
the  integral  only  once  for  each  region. 

Since  the  area  of  the  surface  regions  is  proportional  to  r0 
(see  equation  (8)),  with  this  approximation  the  computa¬ 
tional  time  for  the  evaluation  of  the  interaction  between 
two  surfaces  is  inversely  proportional  to  their  relative 
distance  tq.  Moreover,  the  probability  that  one  build¬ 
ing  A\  “can  see  another  building”  A  2  (i-e.  there  aren’t 
buildings  between  A\  and  A2)  decreases  exponentially 
with  vq  .  Then,  it  can  be  demonstrated  that  the  compu¬ 
tational  time  grows  as  AT§  with  the  number  N  of  trans¬ 
parent  buildings  and  linearly  with  the  number  of  opaque 
ones. 

In  order  to  check  these  results,  we  have  performed 
some  computations  on  configurations  of  n2  buildings 
with  n  =  2, . . . ,  6  illuminated  by  the  field  emitted  from 
an  antenna.  The  buildings  have  been  uniformly  distrib¬ 
uted  over  a  square  region  ( n  per  direction)  and  the  an¬ 
tenna  has  been  located  at  the  center  of  the  layout.  The 
results  are  shown  in  figure  6  where  N  =  n2  represents 
the  number  of  buildings  and  the  relative  computational 
time  is  defined  as  the  ratio  Tj^/T^  between  the  compu¬ 
tational  time  required  for  N  buildings  and  that  required 
for  4  buildings.  We  can  see  that  the  scaling  of  time  is 
similar  to  N?  for  a  small  number  of  buildings  and  to  N 
further. 

4.2  Reflections  on  the  ground 

In  the  simulation  of  the  radioelectric  propagation  in  ur¬ 
ban  areas,  a  numerical  treatment  of  the  ground  using 
the  Kirchhoff  integral  is,  from  a  computational  point  of 
view,  expensive. 

We  assume  therefore  that  the  reflection  from  the 
ground  can  be  computed  by  Geometrical  Optics.  This 


N 

Figure  6:  Scaling  of  the  relative  time  Tj\r /T4  for  a  regular 
distribution  of  y/N  X  y/N  buildings.  Black  diamonds  corre¬ 
spond  to  times  measured  for  the  Kirchhoff  integral  method 


assumption  cannot  be  done  for  reflection  from  the  build 
ings5  walls.  We  refer  to  figure  7  for  a  brief  explanation. 
The  effects  of  diffraction,  which  are  neglected  by  GO, 
are  noticeable  only  near  the  light-shadow  boundaries  of 
the  incident  and  reflected  field.  These  boundaries  exist 
only  when  we  deal  with  reflections  from  the  wall.  For 
ground  reflections,  the  incident  and  reflected  fields  have 
no  discontinuities. 

For  the  ground- wall  double  reflection,  the  discontinu¬ 
ity  is  inside  the  domain  (see  figure  8),  but  is  partially 
compensated  by  the  discontinuity  due  to  the  wall-ground 
double  reflection.  Since  the  wall  and  the  ground  form  a 
90  degree  angle,  the  two  discontinuities  coincide,  and  it 
can  be  seen  that  the  values  of  the  two  contributions  are 
similar. 

We  notice  that  in  the  case  of  perfectly  conducting 
walls  and  ground,  the  effects  of  diffraction  are  nonex¬ 
istent.  In  fact  the  two  surfaces  give  only  three  dis¬ 
tinct  contributions,  two  corresponding  to  single  reflec¬ 
tions  on  the  walls,  one  corresponding  to  a  double  reflec¬ 
tion.  Each  contribution  can  be  attributed  to  a  different 
image  source  and  is  accounted  for  by  GO. 

In  figure  9,  points  A  and  B  represent  a  source  and 
a  receiver  respectively.  The  source  can  be  either  the 
antenna  or  any  point  of  a  scattering  surface  while  the 
receiver  can  be  either  any  point  of  a  different  scattering 
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Figure  7:  Light-shadow  boundaries  of  the  incident  and  re¬ 
flected  field 


Figure  8:  Geometrical  description  of  the  ground-wall  dou¬ 
ble  reflection  paths 

surface  or  any  point  where  we  want  to  compute  the  field 
value. 


Figure  9:  Geometrical  scheme  to  illustrate  ground  reflec¬ 
tions  treatment 

For  each  direct  path  r  from  A  to  B,  there  is  another 
path  1*1  and  r2  which  involves  one  reflection  on  the 
ground  at  the  point  C.  The  contribution  of  the  reflec¬ 
tion  in  C  to  the  field  value  in  B  is  computed  considering 
the  lengths  r1?  r2  and  by  using  Fresnel  coefficients. 

These  considerations  have  been  validated  by  imple¬ 
menting  both  the  explicit  (evaluation  of  the  Kirchhoff 
integral)  and  implicit  (Geometrical  Optics)  treatments 


of  reflections  on  the  ground  and  comparing  the  results 
on  a  test  configuration.  The  layout  (position  of  buildings 
and  observation  axis)  is  schematically  described  in  figure 
2.  The  input  data  of  the  simulation  are  equal  to  those 
given  in  section  3,  with  e  =  4.  The  implicit  treatment  of 
ground  reflections  brings  to  numerical  results  which  are 
similar  to  those  obtained  with  an  explicit  one  (see  figure 
10)  and  it  allows,  at  the  same  time,  a  sensible  reduc¬ 
tion  in  both  the  computational  time  and  storage  mem¬ 
ory.  Moreover,  it  avoids  the  generation  at  the  ground 
boundary,  of  fictitious  diffraction  terms  which  decrease 
the  computed  field  when  ground  reflections  are  explicitly 
considered. 


Figure  10:  Component  Ex  of  the  electromagnetic  field. 
The  continuous  line  represents  the  field  obtained  by  explic¬ 
itly  computing  ground  reflections  with  the  Kirchhoff  inte¬ 
gral;  the  dashed  line  represents  the  field  obtained  by  implicit 
(  purely  GO  )  ground  treatment 


5  Results 

In  sections  3  and  4,  we  have  presented  some  results 
related  to  the  convergence  and  time  complexity  of  the 
method  described  here.  The  validation  of  the  approach 
described  in  this  paper  is  done  by  comparing  the  simu¬ 
lation  outcome  with  a  number  of  experimental  measure¬ 
ments.  In  this  section,  we  are  going  to  present  the  results 
of  these  comparisons. 


5.1  Radioelectric  coverage  in  an  urban 
environment 

We  give  here  the  results  of  a  numerical  simulation  of  the 
propagation  of  the  field  emitted  by  an  antenna  through¬ 
out  an  urban  environment.  In  figure  11,  the  white  blocks 
represent  the  buildings  of  a  real  city  and  their  parameters 
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(dielectric  constant  e,  height  /i&,  etc.)  have  been  chosen 
uniformly  over  the  layout  because  of  the  lack  of  detailed 
measurements.  In  particular,  we  have  chosen  e  =  4  and 
hb  =  10  meters.  The  frequency  of  the  field  emitted  by 
the  source  (antenna  with  a  nonisotropic  radiation  pat¬ 
tern)  is  1.8  GHz  which  corresponds  to  a  wavelength  of 
16.7  centimeters.  Input  power  is  250  mW.  The  source  is 
located  at  a  height  of  4.5  meters  and  the  electric  field  is 
computed  at  a  height  of  1.6  meters. 

The  cross  in  the  center  of  the  figure  marks  the  position 
of  the  source  and  the  gray  scale  indicates  the  simulated 
field  power,  decreasing  from  white  to  black.  For  an  area 
of  about  300  meters  diameter,  the  whole  simulation  has 
taken  24  hours  on  a  RISC/6000  560  IBM. 


Figure  11:  Example  of  the  power  emitted  by  an  anisotropic 
antenna,  obtained  by  the  Kirchhoff  integral  method. 


5.2  Comparison  with  experimental  mea¬ 
surements 

Comparisons  with  experimental  measurements  in  realis¬ 
tic  conditions  are  quite  difficult  because  of  the  absence  of 
accurate  data  on  the  computational  domain.  The  dielec¬ 
tric  constant  of  the  involved  materials  and  the  height  of 
the  buildings  are  not  known.  Moreover,  the  map  of  the 
town  is  available  with  a  precision  of  one  meter;  therefore, 
the  direction  error  for  shorter  walls  can  be  of  the  order 
of  about  ten  degrees  which  doubles  when  we  consider  the 
direction  of  the  reflected  field.  Further,  due  to  interfer¬ 
ence  phenomena,  the  field  can  vary  considerably  within 
one  meter.  Finally,  no  information  is  available  about 


the  presence  of  objects  other  than  buildings;  a  parked 
car  can  easily  intercept  the  ground-reflected  field. 

In  these  conditions,  the  computed  power  may  easily 
differ  by  10  to  20  dbm  from  the  measured  one  and  point- 
to-point  comparisons  cannot  be  used  to  validate  or  in¬ 
validate  the  code.  On  the  other  hand,  the  errors  are 
fairly  localized  and  the  code  can  be  still  used  with  some 
confidence  to  estimate  the  coverage  areas. 

In  the  following  figures,  the  cross  marks  the  position 
of  the  antenna,  while  the  numbered  dots  indicate  the 
positions  of  the  points  where  measurements  have  been 
taken.  The  light  gray  defines  the  area  where  the  com¬ 
puted  power  is  above  threshold  for  both  Geometrical  Op¬ 
tics  and  Kirchhoff’s  method.  The  dark  gray  indicates  the 
area  where  the  power  values  are  above  threshold  only 
for  Kirchhoff’s  method.  The  area  where  the  computed 
power  values  are  below  threshold  is  represented  in  black. 
The  input  data  for  both  simulations  (buildings’  parame¬ 
ters,  source’s  frequency,  area’s  diameter)  are  equal  to 
those  used  in  section  5.1.  Experimental  data  were  ob¬ 
tained  by  a  DECT  system  using  propagation  tester  Sym- 
bionics  SP935. 

In' figure  12,  the  computed  coverage  area  of  the  an¬ 
tenna  is  represented  for  a  threshold  power  of -85  dbm.  It 
can  be  seen  that  the  coverage  ^area  evaluated  with  Geo¬ 
metrical  Optics  is  about  20%  less  than  obtained  with 
Kirchhoff’s  method.  Moreover,  two  of  the  points  consid¬ 
ered  ({5,7})  lie  outside  it,  while  all  the  measured  values 
and  those  evaluated  with  Kirchhoff’s  method  are  above 
threshold.  Figure  13  has  been  obtained  by  raising  the 
threshold  up  to  -70  dbm.  The  measured  values  fall  under 
the  threshold  in  four  points  {1, 2,  3, 4}  and  three  of  them 
are  outside  or  very  close  to  the  boundaries  of  Kirchhof¬ 
f’s  coverage  area.  The  only  evident  discrepancy  concerns 
point  {4}  whose  measured  value  is  under  the  threshold 
while  both  Geometrical  Optics  and  Kirchhoff  s  calcula¬ 
tions  locate  it  inside  the  illuminated  region.  Therefore, 
on  the  one  hand,  Kirchhoff’s  estimate  of  the  coverage 
area  clearly  improves  the  Geometrical  Optics  one;  on 
the  other  hand,  due  to  the  inaccuracy  of  the  input  data, 
the  reliability  of  the  results  is  limited  and  it  is  advisable 
to  use  a  fairly  large  margin  of  safety. 


6  Conclusions 

The  starting  point  for  the  planning  and  maintenance  of 
any  mobile  communication  network  is  based  on  a  pre¬ 
diction  of  the  radio  wave  attenuation  around  the  source 
(antenna)  in  the  layout  of  the  buildings.  In  this  paper 
we  have  proposed  an  iterative  procedure  based  on  the 
Kirchhoff  integral  method.  Results  have  shown  that  this 
method  is  well  suited  to  solve  this  problem  and  offers 
promising  prospects  of  development  in  this  field.  The 
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Figure  12:  The  coverage  region  of  the  antenna  fora  thresh¬ 
old  power  of -85  dbm.  The  light  gray  area  defines  the  zone 
where  both  Geometrical  Optics  and  Kirchhoffs  method 
predict  a  value  for  the  field  power  above  threshold.  The 
dark  gray  area  corresponds  to  the  zone  where  only  the  value 
predicted  by  Kirchhoff’s  method  is  above  threshold.  The 
black  area  defines  the  zone  where  both  Geometrical  Op¬ 
tics  and  Kirchhoff’s  method  predict  values  below  threshold. 
The  antenna  location  is  indicated  by  the  cross,  while  the 
numbered  dots  indicate  the  location  of  receivers  for  exper¬ 
imental  measurements 

model  is  currently  under  revision  to  improve  its  time 
requirements  (Pisani,  1997). 

7  Acknowledgements 

We  would  like  to  thank  Prof.  G.  Mazzarella  of  the  Uni¬ 
versity  of  Cagliari  for  helpful  discussions  and  comments 
on  the  manuscript.  We  acknowledge  A.  Loche  for  provid¬ 
ing  us  with  a  graphic  interface  for  the  visualization  of  the 
radioelectric  coverage  and  Ing.  F.  Tallone  of  the  Centro 
Studi  e  Laboratori  Telecomunicazioni  (CSELT)  of  Turin 
for  providing  us  with  experimental  measurements.  This 
work  has  been  carried  out  with  the  financial  support  of 
Telecom  Italia  and  the  Sardinia  Region  Authorities. 

References 

•  F.  Ikegami,  T.  Takeuchi,  and  S.  Yoshida,  “Theoret¬ 
ical  prediction  of  mean  field  strength  for  urban  mo¬ 
bile  radio”,  IEEE  Trans.  Antenn.  Propagat., 


Figure  13:  The  coverage  region  of  the  antenna  for  a  thresh¬ 
old  power  of  -70  dbm.  The  meaning  of  the  symbols  and 
grey  scale  is  the  same  as  figure  12 

vol.  39,  pp.  299-302,  1991. 

•  J.  D.  Jackson,  Classical  Electrodynamics,  John 
Wiley  &  Sons,  New  York  (1975). 

•  T.  Kiirner,  D.J.  Cichon,  and  W.  Wiesbeck,  “Con¬ 
cepts  and  results  for  3D  digital  terrain-based  wave 
propagation  models:  an  overview”,  IEEE  J.  Se¬ 
lect.  Areas  Commun.,  vol.  11,  no.  7,  pp.  1002- 
12,  1993. 

•  P.  0.  Luthi,  B.  Chopard,  J-F  Wagen,  “  Wave  Prop¬ 
agation  in  Urban  Microcells:  a  Massively  Paral¬ 
lel  Approach  Using  the  TML  Method”,  Lecture 
Notes  in  Computer  Science,  vol.  1041,  pp.  408, 
1996. 

•  L.  Pisani,  “  Dimensionality  reduction  in  the  numer¬ 
ical  evaluation  of  electro-magnetic  scattering:  from 
the  Kirchhoff’s  integral  method  to  a  zero  dimen¬ 
sional  approach”,  submitted  to  Computer  Phys. 
Commun.. 

•  K.  Rizk,  J-F  Wagen,  F.  Gardiol,  “  Ray  tracing  based 
path  loss  prediction  in  two  microcellular  environ¬ 
ments”,  in  Proceedings  IEEE  PIMRC’94,  pp. 
210-4,  The  Hague,  Netherlands,  1994. 

•  J.P.  Rossi,  A.J.  Levy,  “A  ray  model  for  decimetric 
radiowave  propagation  in  an  urban  area”,  Radio 
Science,  vol.  27,  no.  6,  pp.  971-9,  1992. 


70 


On  the  Bounded  Part  of  the  Kernel  in  the 
Cylindrical  Antenna  Integral  Equation 

M.P.  Ramachandran 

Analysis  and  Facility  Division 
Space  Applications  Center 
Ahmedabad  380  053,  India 


Abstract 

The  kernel  in  the  cylindrical  antenna  integral  equation  was  partitioned  by  Schelkunoff  into  a  com¬ 
plete  elhptic  integral  and  a  bounded  integral.  This  paper  gives  an  exact  expression  for  the  bounded 
part. 

1  Introduction 

The  cylindrical  antenna  integral  equation  for  total  axial  current  I (z)  on  perfectly  conducting  tube 
of  length  ‘2h’  and  radius  ‘a’  is  [1]: 

ju(J^  +  k?jA,(z)  =  k'2Ei(z),  (1) 

where  h 

Mz)  =  £  j_h  K(Z  -  z')I(z')dz' .  (2) 

The  kernel  in  (2)  is: 

K{z-£)=h 

where  x,2 

R  =  R(z  —  z',  <f>')  =  (z  -  z'f  +  4a2  sin2  (yl  .  (4) 

Here  e  and  p.  characterize  the  medium  [1],  and  k  is  the  wavenumber. 

For  thin  structures  (a  «  h  and  a  «  A),  K(z  -  z')  is  usually  substituted  by  a  reduced  kernel 
approximation  ^ 

K~  Kr{z-z')=(—— ,  (5) 

r 

where  , 

r  =  r(z  —  z')  —  (z  —  z'f  +  a2j  .  (6) 

The  kernel  (3)  possesses  a  logarithmic  singularity  (see  for  example  [2]).  This  needs  to  be  included 
while  solving  (1).  Schelkunoff  has  partitioned  the  kernel  into  two  parts. 

The  first  part  is  a  complete  elliptic  integral  and  the  second  part  is  a  bounded  integral.  By 
expanding  the  elliptic  integral,  Pearson  [3]  has  derived  an  expression  which  includes  the  logarithmic 
singularity.  On  the  other  hand,  as  the  exact  expression  of  the  bounded  integral  is  not  available, 
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it  is  usually  evaluated  numerically.  Karwowski  [4]  has  suggested  a  closed  form  approximation;  an 
alternative  to  the  numerical  integration.  In  this  presentation  we  obtain  an  exact  expression  for  the 
bounded  part.  Certain  recurrent  relations  are  obtained  to  simplify  its  computation.  The  stability 
of  the  recurrence  relations  is  analyzed  and  checked. 

Before  we  conclude  this  section  we  shall  briefly  review  other  research  accomplishments  regarding 
the  kernel  (3).  Wang  [5]  has  obtained  an  ‘exact’  expression  for  the  kernel  (3)  in  terms  of  spherical 
Hankel  functions.  Werner  [6]  has  presented  an  alternative  ‘exact’  expression  replacing  the  spherical 
Hankel  functions  by  complex  exponential  functions.  This  has  an  advantage  from  the  point  of  view  of 
analytical  and  numerical  evaluation  of  the  kernel.  Werner  has  presented  ‘extended’  approximations 
[7]  to  the  kernel,  by  truncating  the  series  to  a  few  terms.  This  is  an  alternative  to  the  reduced 
approximation  (5).  Interestingly,  the  vector  potential  in  (2)  due  to  the  singular  part  KB  has  been 
worked  by  Werner,  et  al  [8].  The  results  of  this  paper  on  the  bounded  part  Kb  can  be  useful  when 
combined  with  the  results  of  Pearson  [3]  or  Werner,  et  al  [8]. 


2  Kernel 


Recall  the  Schelkunoff  partition: 
where 

and 


K{z-z')  =  Ke  +  Kb, 

1  rn  I 

*•'2 ~,Lr^ 

1  n-e-*V 


**  =  -2 -J- 


R 


(V 

(8) 

(9) 


Expanding  the  first  integral,  an  elliptic  integral,  Pearson  has  obtained  the  following  expression, 
which  we  reproduce,  with  a  correction: 


where 


K, 


EB 


Ke  = - In  \z  —  z'\  +  Keb  , 

Tva 


1-/?,  ,  P 

- In \z  —  2d - - 

7 ra  7r  a 


(10) 


'IN  2 

In  8a  —  In  /3  +(  - 

\  Zj  , 


In 


+(r-D[ta(s)-(r5)-(r-J 

We  have  used  the  following  notation 


A. 

2 


-(rs) 


Pi 

Pi- 


p  = 


2  a 


[4a2  +  (z-  z')2]l/2 


,  A  =  (1  -  p2)1/2 . 


(11) 


(12) 


The  series  is  convergent  for  \z  -  z'\  <  4a.  The  leading  term  in  KEB  can  be  shown  to  be  of  the  form 
(z  —  z')2  In  | z  —  z'\.  Next,  we  shall  obtain  an  exact  expression  for  KB. 
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(13) 


3  Expression  for  Kb 

Kb  can  be  written  as 

Kb  =  Kbt  +  &Bi  • 

where  KBt  is  the  real  part  and  KBi  is  the  imaginary  part 


—1  /-7r  (1  -  cos  kR)  ,,, 

K°'=2*L  R  « 

(14) 

and 

-j  r*  sinfci?,., 

■**“» *L  R  d *■ 

(15) 

Expanding  sin  kR  in  eq.  (15)  and  setting  ( z  —  zr )  —  u,  we  get  the  infinite  series. 

j  ^  (—l)mk2m+1  Am-i(u,  a) 

Bi~  (2m  +  1)! 

(16) 

where 

Am_i(u,a)  =  /  Rlmdfx 

Jo 

(17) 

and 

Ri  =  [u2  +  4a2  sin2  /ij  /  . 

(18) 

We  find  from  eqs.  (17)  and  (18)  that 

A_i(w,  a)  =  7r 

(19) 

and 

a)  =  7r  fu2  +  2a2")  . 

(20) 

In  eq.  (17),  we  rewrite  Am-i(u,a )  as: 

Am-\ (u,  a)  =  r  (u2  +  4 a2  sin2  ju)™  1  (a2  +  4a2  sin2  //)  d\x 

=  (a2  +  2a2)  Am_2(u,  a)  -  2 a2  /  (a2  +  4a2  sin2  /i) ”*  1  cos  2/ui/z . 

Apply  integration  by  parts  to  the  integral.  After  certain  algebraic  manipulations,  we  get  the  recur¬ 
rence  relations: 

[1  +  (m  —  1)]  Am_ i  =  (u2  +  2a2)  [1  +  2 (m  —  1)]  Am_2 

-(m  -  l)u2  (a2  +  4a2)  Am_3  ;  m  =  2, 3, 4,  -  -  •  (21) 

Using  A_i,  A0  from  eqs.  (19)  and  (20)  we  find  Am-\\  m  =  2, 3, 4,  •  •  •  and  then  determine  KBi  from 
eq.  (16)  up  to  a  desired  accuracy.  The  stability  of  the  recurrence  relations  needs  to  be  examined 
for  a  given  radius  ‘a’  and  length  ‘2h’  of  the  wire.  This  has  been  done  in  the  stability  section.  Now, 
on  similar  lines  the  real  part,  namely  Kqt ,  becomes 

KB'  =  ~Ik  (2m  +  2)!  (22) 
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where 


Bm-i(u,a)=  I  Rjm+1dfi. 
Jo 


We  discuss  the  computation  of  Kbt  for  two  distinct  cases  of  u  =  z  —  z',  namely  for  u  =  0  and 
for  m^O. 

Prom  eq.  (23)  we  find 

B_i(0,o)  =  4a  (24) 

and  then  obtain  recurrence  relation 


Bm- i(0,o)—  (27ft  +  1)  Bm— 2(0,0).  (25) 

Computed  values  of  Bm-i,  m  =  1,2,3,  ••  •  are  substituted  in  eq.  (22).  It  may  be  noted  that  the 
recurrence  relation  is  stable  for  a  <  0.5.  This  is  a  dimensionless  number. 

We  obtain  the  following  results  using  [9] 


/*7r/2  ,  .  1/2 

5_i(u,  o)  =  2  uj  (l  +  g2sin/ij  dfi 


2 usE  ( ^ 


where  we  have  used  the  notation  q  =  (2 a/u),  s  —  (1  +  q2)1/2.  Further, 


B0(u,  a)  =  2 v?  \sE  f  +  q2I(q)  , 


where 


E  and  F  are  elliptic  integrals  of  second  and  first  kind,  respectively,  and  are  evaluated  to  a  desired 
accuracy  using  Landen’s  Transformation  [10].  The  recurrence  relation  now  is: 


Bra — 1  (u,  a)  1  + 


2m  —  1 


—  (u2  +  2 a2^  Bm- 2(w,  a)  [1  +  (2m  —  1)] 


2m  —  1 


u2  (u2  +  4 a2)  5m_ 3(71,  a) 


m  =  2, 3, 4,  •  •  •  (29) 

Using  eqs.  (26)  and  (27)  we  can  obtain  Bm-  1,  m  =  2, 3, 4,  •  •  •  and  then  evaluate  The  stability 
of  eq.  (29)  is  also  included  in  the  next  section. 


4  Stability 

The  characteristic  equation  [11]  corresponding  to  eq.  (21)  is: 

(1  +  n)x 2  -  <*(1  +  2n):r  +  717  =  0  ;  n  =  1, 2,  •  •  • 
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(31) 


where  a  =  (u2  +  2a2),  7  =  u2  ( u 2  +  4a2).  Rewrite  eq.  (30)  as 


x(x  —  a)  _  n 
ax  —  7  (1  +  n) 


The  above  equation  is  useful  to  seek  the  stability  condition,  namely  \x\  <  1.  For  a  given  radius  ‘a’ 
and  length  of  wire  ‘2 h\  this  condition  shall  depend  upon  A. 


Denoting 


n 

(l  +  n) 


(32) 


The  solution  of  (31)  is 

_  a(l  +  P)  ±  [a2(l  ~  Pf  +  16pa4]1/2  (33\ 

x~  2  '  K 
The  right  hand  side  is  a  monotonic  increasing  function  for  given  values  of  p  and  a,  enabling  direct 
evaluation  of  the  maximum  of  x.  For  recurrence  relations  corresponding  to  the  real  part  given  in 
eq.  (29),  the  analysis  is  repeated  on  similar  lines  with 


(2  n  +  1) 

P  =  (2n  +  2)  ’ 


n  —  1, 2, 3,  *  •  *  , 


(34) 


yet  the  stability  criteria  is  not  altered.  Thus,  for  any  wire  of  known  radius  ‘a’  and  length  ‘2 h\  using 
eq.  (33)  we  can  easily  establish  the  condition  of  stability  on  A  before  computing  the  recurrence 
relations  and  determining  Kb- 


5  Examples 

As  a  first  case  we  have  a  dipole  of  half-wave  length  in  length  and  whose  radius  is  (A/8).  From  eq. 
(21),  we  note  that  the  value  of  p  varies  from  1/2  to  a  limit  of  1.  Then,  the  maximum  of  x  in  eq.  (33) 
is  5A2/16,  which  implies  that  the  condition  for  stability  is  A  <  4/\/5,  that  is  approximately  1.7888. 
This  is  not  a  strong  condition  as  only  a  limiting  value  could  be  used.  In  the  first  example  we  set  the 
radius  of  the  wire  to  be  0.003  and  A  to  be  1.0.  This  is  a  case  of  thin  wire  and  Karwowski’s  results 
are  accurate.  We  have  given  in  Table  1  for  various  values  of  u  =  z  —  z' ,  the  (real,  imaginary)  values 
of  Kb  from  (16)  and  (22).  These  values  are  compared  with  those  of  Karwowski  [4], 

Next,  we  consider  a  cylindrical  wire  of  a  wavelength  in  length  and  radius  (A/4).  The  maximum 
value  of  x  is  (5A2/4).  Thus,  for  A  <  2/\/5  (approximately  0.8944),  the  recurrence  relations  in  eqs. 
(21)  and  (29)  are  stable.  The  condition  in  eq.  (25)  that  a  <  1/2  is  satisfied  in  this  case  also.  In 
the  second  example  we  have  a  =  0.22  and  A  =  0.88.  The  corresponding  results  are  given  in  Table 
2.  Here,  while  setting  A  =  0.90,  the  recurrence  relation  was  verified  to  be  unstable,  which  agrees 
with  our  estimate.  The  values  in  the  tables  are  correct  up  to  four  decimal  places  and  appropriately 
rounded  up  to  three  decimal  places.  Only  12  terms  in  the  series  were  used  to  obtain  this  accuracy. 
The  difference  between  the  sets  of  the  values  is  due  to  the  approximation  introduced  by  Karwowski 
in  replacing  (1  / R)  by  (sin  4>/R)  in  (9)  and  thus  apparently  diminishing  the  contribution  near  the 
ends  of  the  interval. 
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Table  1:  Comparison  between  this  work  and  that  of  Karwowski  [4]  for  a  thin  wire  half-wave  dipole. 


a  = 

0.003  and 

A  =  1.0 

u 

Present 

Karwowski 

0.0 

(-0.075, 

-6.282) 

(-0.075, 

-6.282) 

0.1 

(-1.885, 

-5.877) 

(-1.911, 

-5.877) 

0.2 

(-3.433, 

-4.755) 

(-3.455, 

-4.755) 

0.3 

(-4.347, 

-3.170) 

(-4.364, 

-3.170) 

0.4 

(-4.513, 

-1.469) 

(-4.552, 

-1.469) 

0.5 

(-3.997, 

+0.226) 

(-4.000, 

+0.226) 

Table  2:  Comparison  between  this  work  and  that  of  Karwowski  [4]  for  a  thick  wavelength  long  wire. 


u 

a  =  0.22  and  A 

Present 

=  0.88 

Karwowski 

0.000 

(-4.114,  -3.063) 

(-4.341,  -2.894) 

0.176 

(-4.500,  -2.001) 

(-4.584,  -1.852) 

0.352 

(-4.150,  +0.150) 

(-4.152,  +0.250) 

0.528 

(-2.577,  +1.406) 

(-2.215,  +1.451) 

0.704 

(-0.884,  +0.903) 

(-0.431,  +0.909) 

0.880 

(-0.323,  -0.371) 

(-0.096,  -0.380) 

6  Conclusion 

This  paper  makes  an  effort  to  evaluate  the  kernel  by  the  definition  of  Schelkunoff  [1].  An  exact 
expression  for  the  bounded  part  of  the  kernel  KB  is  given.  This  result  can  be  used  along  with  that 
of  Pearson  [3]  for  solving  cylindrical  antenna  integral  equations.  Examples  are  considered  where 
Kb  is  evaluated  by  the  present  method  and  is  compared  with  that  of  Karwowski  [4] . 
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LETTER  FROM  THE  EDITOR 


Dear  colleagues  and  friends  in  the  ACES  Community 

Sadly,  life  moves  on  and  there  comes  a  time  when  we  shed  some  responsibilities  and  take  on  new 
challenges.  This  will  be  my  last  issue  of  the  ACES  Journal  as  Editor-in-Chief.  Our  very  able  colleagues 
Allen  Glisson  and  Ahmed  Kishk  will  be  at  the  helm  of  the  ACES  Journal  for  the  next  few  years. 

When  I  initially  accepted  the  responsibility  for  the  Journal,  it  was  with  some  trepidation  because  of 
possible  distance  and  communications  problems.  On  the  whole  these  fears  have  proved  groundless  because 
of  the  extensive  use  made  of  email  and  fax.  Yes,  there  have  been  frustrations.  Murphy’s  law  dictates  that 
our  email  servers  must  crash  whenever  a  printing  deadline  is  in  the  offing.  Believe  me,  you  really  do  not 
want  to  sit  with  a  15  January  deadline  over  the  peak  holiday  season  in  South  Africa  when  everything  pretty 
well  grinds  to  a  halt.  Incredible  as  it  may  seem,  of  all  the  hundreds  of  pieces  of  mail,  only  a  very  few 
disappeared  forever  en  route  to  reviewers  (fortunately  never  to  authors).  Airmail  delivery  from  South  Africa 
to  North  and  South  America,  Europe  and  other  distant  parts  has  also  been  known  to  take  up  to  6  weeks!  I 
have  no  idea  where  the  problems  lie. 

Handing  over  editorial  responsibilities  is  a  drawn-out  process.  Many  of  the  papers  currently  in  the 
review  cycle  will  be  published  under  Allen’s  and  Ahmed’s  stewardship.  Hopefully  they  will  benefit  as  much 
from  my  efforts  as  I  did  from  my  predecessor  David  Stein  when  I  started  taking  over  editorship  for  the 
Journal  in  July/August  of  1993.  At  the  outset  I  had  the  wholehearted  support  and  guidance  of  the  ACES 
Editor-in-Chief,  Perry  Wheless,  and  also  the  able  assistance  and  support  of  Bela  Konrad,  the  Journal’s 
Associate  Editor-in-Chief.  It  has  been  a  particular  privilege  to  be  able  to  work  so  closely  with  them  over  the 
intervening  years. 

Inevitably  a  new  broom  wants  to  bring  in  changes  and  improvements.  The  ACES  Committee 
supported  the  decision  to  move  to  the  double  column  compact  format.  This  has  resulted  in  significant 
savings  in  publications,  and  allowed  us  more  flexibility  as  regards  limitations  on  paper  lengths.  However, 
the  one  thing  I  had  hoped  to  achieve  has  eluded  me.  Hopefully  my  successors  will,  with  your  support, 
succeed  in  making  the  Journal  a  quarterly  publication. 

Collating  the  efforts  of  authors,  editors  and  reviewers  would  however  be  a  futile  exercise  without 
someone  to  handle  the  final  publication  and  distribution,  not  only  of  the  Journal  but  also  the  Newsletter. 
Here  we  in  ACES  are  particularly  fortunate  to  have  the  services  of  the  remarkable  husband  and  wife  team, 
Dick  and  Pat  Adler.  Dick  of  course  has  been  Managing  Editor  for  ACES  since  its  earliest  days.  Not  many  of 
our  members  are  aware  of  the  behind  the  scenes  activities  and  contributions  of  these  two  -  and  that  is  as  it 
should  be  in  as  well  rim  and  smooth  an  operation  as  theirs.  The  reader  gets  a  quality  product  because  Dick 
and  Pat  set  a  high  standard.  On  more  than  one  occasion  this  has  required  redoing  paste-ups  for  articles, 
reprinting  the  master  copy  from  one  of  many  incompatible  wordprocessing  formats,  and  more.  Pat  and 
Dick,  on  behalf  of  our  members  and  myself,  “Many  thanks.”  Without  your  support  I  would  not  have  coped. 

At  the  end  of  this  letter,  I  will  again  be  thanking  many  of  our  colleagues  who  are  not  on  the  editorial 
panel  shown  in  the  Journal.  I  acknowledge  their  willing  support  and  assistance  over  the  years  with  a  great 
deal  of  gratitude.  The  bulk  of  the  reviews  were  undertaken  by  our  international  panel.  I  thank  them  for 
their  assistance,  their  patience  and  forbearance,  and  particularly  their  time.  All  out  reviewers  and  editors 
have  heavy  professional  commitments.  That  they  do  find  and  make  time  for  the  reviews  is  a  tribute  to  their 
commitment  to  ACES  and  its  members. 

The  Guest  Editors  of  the  various  Special  Issues  are  indeed  a  Special  breed.  They  have  taken  on  an 
additional  volunteer  activity  on  top  of  other  commitments.  For  their  dedication  and  many  hours  of  hard 
work  and  excellent  contributions,  I  thank  them  on  behalf  of  all  ACES  members. 
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As  I  have  said  on  more  than  one  occasion,  the  Journal  really  belongs  to  our  members,  and  it  is  the 
authors  who  make  it  possible.  I  have  no  doubt  that  you  will  give  the  new  editors  your  wholehearted 
support,  and  I  look  forward  to  reading  many  more  fascinating  articles  covering  a  wide  diversity  of  interests 
and  applications  in  the  field  of  computational  electromagnetics. 

Lastly,  it  is  my  pleasure  to  acknowledge  the  support  of  the  Department  of  Electrical  and  Electronic 
Engineering  at  the  University  of  Pretoria.  Not  only  has  it  provided  email,  mail  and  fax  facilities,  but  it  has 
also  made  available  the  assistance  of  two  of  the  administrative  staff,  Ms.  Alet  van  Zyl,  and  Ms.  Cornel 
Freislich.  Alet  and  Cornel,  my  sincere  thanks  to  you  both. 

Editorship  of  the  ACES  Journal  has  been  a  stimulating  and  productive  part  of  my  professional 
career.  Thank  you  for  the  opportunity  to  be  of  service  to  you,  the  ACES  Member.  It’s  been  fun. 

Duncan  C.  Baker 

Outgoing  Editor-in-Chief,  ACES  Journal 
16  January  1998. 
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disclosure  of  their  material  requires  the  prior  consent  of  other  parties  and  if  so,  to  obtain  it.  Furthermore,  ACES  must  assume  that,  if  an  author 
uses  within  his/her  article  previously  published  and/or  copyrighted  material  that  permission  has  been  obtained  for  such  use  and  that  any  required 
credit  lines,  copyright  notices,  etc.  are  duly  noted. 

AUTHOR/COMPANY  RIGHTS 

If  you  are  employed  and  you  prepared  your  paper  as  a  part  of  your  job,  the  rights  to  your  paper  initially  rest  with  your  employer.  In  that  case, 
when  you  sign  the  copyright  form,  we  assume  you  are  authorized  to  do  so  by  your  employer  and  that  your  employer  has  consented  to  all  of  the 
terms  and  conditions  of  this  form.  If  not,  it  should  be  signed  by  someone  so  authorized. 

NOTE  RE  RETURNED  RIGHTS:  Just  as  ACES  now  requires  a  signed  copyright  transfer  form  in  order  to  do  “business  as  usual”, 
it  is  the  intent  of  this  form  to  return  rights  to  the  author  and  employer  so  that  they  too  may  do  “business  as  usual”.  If  further  clarification  is  required, 
please  contact:  The  Managing  Editor,  R.W.  Adler,  Naval  Postgraduate  School,  Code  EC/AB,  Monterey,  CA,  93943,  USA  (408)656-2352. 

Please  note  that,  although  authors  are  permitted  to  re-use  all  or  portions  of  their  ACES  copyrighted  material  in  other  works,  this  does  not 
include  granting  third  party  requests  for  reprinting,  republishing,  or  other  types  of  re-use. 

JOINT  AUTHORSHIP 

For  jointly  authored  papers,  only  one  signature  is  required,  but  we  assume  all  authors  have  been  advised  and  have  consented  to  the  terms 
of  this  form. 

U.S.  GOVERNMENT  EMPLOYEES 

Authors  who  are  U.S.  Government  employees  are  not  required  to  sign  the  Copyright  Transfer  Form  (Part  A),  but  any  co-authors  outside  the 
Government  are. 

Part  B  of  the  form  is  to  be  used  instead  of  Part  A  only  if  all  authors  are  U.S.  Government  employees  and  prepared  the  paper  as  part  of  their 
job. 

NOTE  RE  GOVERNMENT  CONTRACT  WORK:  Authors  whose  work  was  performed  under  a  U.S  Government  contract  but  who  are 
not  Government  employees  are  required  to  sign  Part  A-Copyright  Transfer  Form.  However,  item  5  of  the  form  returns  reproduction  rights  to 
the  U.S.  Government  when  required,  even  though  ACES  copyright  policy  is  in  effect  with  respect  to  the  reuse  of  material  by  the  general  public. 
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INFORMATION  FOR  AUTHORS 


PUBLICATION  CRITERIA 

Each  paper  is  required  to  manifest  some  relation  to  applied 
computational  electromagnetics.  Papers  may  address 
general  issues  in  applied  computational  electromagnet¬ 
ics,  or  they  may  focus  on  specific  applications,  tech¬ 
niques,  codes,  or  computational  issues.  While  the 
following  list  is  not  exhaustive,  each  paper  will  generally 
relate  to  at  least  one  of  these  areas: 

1.  Code  validation.  This  is  done  using  internal  checks  or 
experimental,  analytical  or  other  computational  data. 
Measured  data  of  potential  utility  to  code  validation  efforts 
will  also  be  considered  for  publication. 

2.  Code  performance  analysis.  This  usually  involves 
identification  of  numerical  accuracy  or  other  limitations, 
solution  convergence,  numerical  and  physical  modeling 
error,  and  parameter  tradeoffs.  However,  it  is  also 
permissible  to  address  issues  such  as  ease-of-use,  set-up 
time,  run  time,  special  outputs,  or  other  special  features. 

3.  Computational  studies  of  basic  physics.  This  involves 
using  a  code,  algorithm,  or  computational  technique  to 
simulate  reality  in  such  a  way  that  better  or  new  physical 
insight  or  understanding  is  achieved. 

4.  New  computational  techniques,  or  new  applications  for 
existing  computational  techniques  or  codes. 

5.  ’’Tricks  of  the  trade"  in  selecting  and  applying  codes 
and  techniques. 

6.  New  codes,  algorithms,  code  enhancement,  and  code 
fixes.  This  category  is  self-explanatory  but  includes 
significant  changes  to  existing  codes,  such  as  applicability 
extensions,  algorithm  optimization,  problem  correction, 
limitation  removal,  or  other  performance  improvement. 
Note:  Code  (or  algorithm)  capability  descriptions  are 
not  acceptable,  unless  they  contain  sufficient  technical 
material  to  justify  consideration. 

7.  Code  input/output  issues.  This  normally  involves 
innovations  in  input  (such  as  input  geometry 
standardization,  automatic  mesh  generation,  or  computer- 
aided  design)  or  in  output  (whether  it  be  tabular,  graphical, 
statistical,  Fourier-transformed,  or  otherwise  signal- 
processed).  Material  dealing  with  input/output  database 
management,  output  interpretation,  or  other  input/output 
issues  will  also  be  considered  for  publication. 

8.  Computer  hardware  issues.  This  is  the  category  for 
analysis  of  hardware  capabilities  and  limitations  in  meeting 


various  types  of  electromagnetics  computational  require¬ 
ments.  Vector  and  parallel  computational  techniques  and 
implementation  are  of  particular  interest. 

Applications  of  interest  include,  but  are  not  limited  to, 
antennas  (and  their  electromagnetic  environments), 
networks,  static  fields,  radar  cross  section,  shielding, 
radiation  hazards,  biological  effects,  electromagnetic  pulse 
(EMP),  electromagnetic  interference  (EMI),  electromagnet¬ 
ic  compatibility  (EMC),  power  transmission,  charge 
transport,  dielectric  and  magnetic  materials,  microwave 
components,  MMIC  technology,  remote  sensing  and  geo¬ 
physics,  communications  systems,  fiber  optics,  plasmas, 
particle  accelerators,  generators  and  motors,  electromagnet¬ 
ic  wave  propagation,  non-destructive  evaluation,  eddy 
currents,  and  inverse  scattering. 

Techniques  of  interest  include  frequency-domain  and 
time-domain  techniques,  integral  equation  and  differential 
equation  techniques,  diffraction  theories,  physical  optics, 
moment  methods,  finite  differences  and  finite  element 
techniques,  modal  expansions,  perturbation  methods,  and 
hybrid  methods.  This  list  is  not  exhaustive. 

A  unique  feature  of  the  Journal  is  the  publication  of 
unsuccessful  efforts  in  applied  computational 
electromagnetics.  Publication  of  such  material  provides  a 
means  to  discuss  problem  areas  in  electromagnetic  model¬ 
ing.  Material  representing  an  unsuccessful  application  or 
negative  results  in  computational  electromagnetics  will  be 
considered  for  publication  only  if  a  reasonable  expectation 
of  success  (and  a  reasonable  effort)  are  reflected. 
Moreover,  such  material  must  represent  a  problem  area  of 
potential  interest  to  the  ACES  membership. 

Where  possible  and  appropriate,  authors  are  required  to 
provide  statements  of  quantitative  accuracy  for  measured 
and/or  computed  data.  This  issue  is  discussed  in  "Accuracy 
&  Publication:  Requiring  quantitative  accuracy  statements 
to  accompany  data",  by  E.K.  Miller,  ACES  Newsletter ,  Vol. 
9,  No.  3,  pp.  23-29,  1994,  ISBN  1056-9170. 

EDITORIAL  REVIEW 

In  order  to  ensure  an  appropriate  level  of  quality  control, 
papers  are  refereed.  They  are  reviewed  both  for  technical 
correctness  and  for  adherence  to  the  listed  guidelines 
regarding  information  content.  Authors  should  submit  the 
initial  manuscript  in  draft  form  so  that  any  suggested 
changes  can  be  made  before  the  photo-ready  copy  is 
prepared  for  publication. 


STYLE  FOR  CAMERA-READY  COPY 

The  ACES  Journal  is  flexible,  within  reason,  in  regard  to 
style.  However,  certain  requirements  are  in  effect: 

1.  The  paper  title  should  NOT  be  placed  on  a  separate 
page.  The  title,  author(s),  abstract,  and  (space  permitting) 
beginning  of  the  paper  itself  should  all  be  on  the  first  page. 
The  title,  author(s),  and  author  affiliations  should  be 
centered  (center-justified)  on  the  first  page. 

2.  An  abstract  is  REQUIRED.  The  abstract  should  state 
the  computer  codes,  computational  techniques,  and 
applications  discussed  in  the  paper  (as  applicable)  and 
should  otherwise  be  usable  by  technical  abstracting  and 
indexing  services. 

3.  Either  British  English  or  American  English  spellings 
may  be  used,  provided  that  each  word  is  spelled  consistently 
throughout  the  paper. 

4.  Any  commonly-accepted  format  for  referencing  is 
permitted,  provided  that  internal  consistency  of  format  is 
maintained.  As  a  guideline  for  authors  who  have  no  other 
preference,  we  recommend  that  references  be  given  by 
author(s)  name  and  year  in  the  body  of  the  paper  (with 
alphabetical  listing  of  all  references  at  the  end  of  the  paper). 
Titles  of  Journals,  monographs,  and  similar  publications 
should  be  in  boldface  or  italic  font  or  should  be  underlined. 
Titles  of  papers  or  articles  should  be  in  quotation  marks. 

5.  Internal  consistency  shall  also  be  maintained  for  other 
elements  of  style,  such  as  equation  numbering.  As  a 
guideline  for  authors  who  have  no  other  preference,  we 
suggest  that  equation  numbers  be  placed  in  parentheses  at 
the  right  column  margin. 

6.  The  intent  and  meaning  of  all  text  must  be  clear.  For 
authors  who  are  NOT  masters  of  the  English  language,  the 
ACES  Editorial  Staff  will  provide  assistance  with  grammar 
(subject  to  clarity  of  intent  and  meaning). 

7.  Unused  space  should  be  minimized.  Sections  and 
subsections  should  not  normally  begin  on  a  new  page. 

MATERIAL,  SUBMITTAL  FORMAT  AND 
PROCEDURE 

The  preferred  format  for  submission  and  subsequent  review, 
is  12  point  font  or  12  cpi,  double  line  spacing  and  single 
column  per  page.  Four  copies  of  all  submissions  should  be 
sent  to  the  Editor-in-Chief  (see  inside  front  cover).  Each 
submission  must  be  accompanied  by  a  covering  letter.  The 
letter  should  include  the  name,  address,  and  telephone 
and/or  fax  number  and/or  e-mail  address  of  at  least  one  of 
the  authors. 


Only  camera-ready  original  copies  are  accepted  for 
publication.  The  term  "camera-ready”  means  that  the 
material  is  neat,  legible,  and  reproducible.  The  preferred 
font  style  is  Times  Roman  10  point  (or  equivalent)  such  as 
that  used  in  this  text.  A  double  column  format  similar  to 
that  used  here  is  preferred.  No  author’s  work  will  be 
turned  down  once  it  has  been  accepted  because  of  an 
inability  to  meet  the  requirements  concerning  fonts  and 
format  Full  details  are  sent  to  the  author(s)  with  the  letter 
of  acceptance. 

There  is  NO  requirement  for  India  ink  or  for  special  paper; 
any  plain  white  paper  may  be  used.  However,  faded  lines  on 
figures  and  white  streaks  along  fold  lines  should  be  avoided. 
Original  figures  -  even  paste-ups  -  are  preferred  over 
"nth-generation"  photocopies.  These  original  figures  will 
be  returned  if  you  so  request. 

While  ACES  reserves  the  right  to  re-type  any  submitted 
material,  this  is  not  generally  done. 

PUBLICATION  CHARGES 

ACES  members  are  allowed  12  pages  per  paper  without 
charge;  non-members  are  allowed  8  pages  per  paper  without 
charge.  Mandatory  page  charges  of  $75  a  page  apply  to  all 
pages  in  excess  of  12  for  members  or  8  for  non-members. 
Voluntary  page  charges  are  requested  for  the  free  (12  or  8) 
pages,  but  are  NOT  mandatory  or  required  for  publication. 
A  priority  courtesy  guideline,  which  favors  members, 
applies  to  paper  backlogs.  Full  details  are  available  from  the 
Editor-in-Chief. 

COPYRIGHTS  AND  RELEASES 

Each  primary  author  must  sign  a  copyright  form  and  obtain 
a  release  from  his/her  organization  vesting  the  copyright 
with  ACES.  Forms  will  be  provided  by  ACES.  Both  the 
author  and  his/her  organization  are  allowed  to  use  the 
copyrighted  material  freely  for  their  own  private  purposes. 

Permission  is  granted  to  quote  short  passages  and  reproduce 
figures  and  tables  from  an  ACES  Journal  issue  provided  the 
source  is  cited.  Copies  of  ACES  Journal  articles  may  be 
made  in  accordance  with  usage  permitted  by  Sections  107 
or  108  of  the  U.S.  Copyright  Law.  This  consent  does  not 
extend  to  other  kinds  of  copying,  such  as  for  general 
distribution,  for  advertising  or  promotional  purposes,  for 
creating  new  collective  works,  or  for  resale.  The 
reproduction  of  multiple  copies  and  the  use  of  articles  or 
extracts  for  commercial  purposes  require  the  consent  of  the 
author  and  specific  permission  from  ACES.  Institutional 
members  are  allowed  to  copy  any  ACES  Journal  issue  for 
their  internal  distribution  only. 


